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I  ABSTRACT 


The  unsteady,  three-dimensional,  incompressible,  viscous  flow  interactions  between  a 
single  vortex  tube  or  a  pair  of  vortex  tubes  advected  by  a  uniform  free  stream  and  a 
spherical  particle  held  fixed  in  space  were  investigated  numerically  for  a  range  of  particle 
Reynolds  numbers  between  20  and  100.  Useful  correlations  of  lift  coefficient,  moment 
coefficient,  and  drag  coefficient  with  velocity  fluctuation,  Reynolds  number,  offset  dis¬ 
tance,  and  initial  vortex  size  were  obtained  and  reported.  A  new  mechanism  based  upon 
droplet  lift  has  been  suggested  for  the  dispersion  of  sprays.  The  mechanism  of  particle 
dispersion  due  to  the  interaction  with  small  vortices  was  quite  different  from  that  due 
to  the  interaction  with  a  large  vortex.  A  new  improved  equation  of  particle  (or  droplet) 
motion  has  been  demonstrated  to  be  superior  to  the  previously  proposed  equations  of 
particle  motion.  The  droplet  heating  study  showed  the  very  significant  response  of  the 
Nusselt  number  to  vortical  disturbances.  The  new  computations  with  thermocapillary 
effects  showed  another  coupling  of  the  fluid  motion  to  the  thermal  field. 


II  OBJECTIVES 

The  objectives  of  this  research  program  axe  to  investigate  the  interactions  of  vapor¬ 
izing  droplets  with  a  turbulent  field  of  the  type  encountered  in  gas  turbine  combustors. 
It  is  intended  to  develop  predictive  capability  through  the  use  of  correlations.  There  is 
special  interest  in  the  important  and  challenging  high-frequency  end  of  the  energy  spec¬ 
trum  where  turbulent  length  scales  are  comparable  to  droplet  size.  The  full  Navier-Stokes 
equations  were  numerically  solved  and  a  simple  mathematical  description  for  the  turbu¬ 
lent  velocity  fluctuation  was  employed.  In  the  mathematical  description,  turbulent-like 
fluctuations  were  simulated  in  a  controlled  way  by  introducing  cylindrical  vortices  which 
have  a  length  scale  of  the  order  of  that  of  the  droplet  and  a  strength  corresponding  to  a 
turbulent  velocity  fluctuation.  From  the  calculations,  instantaneous  lift,  drag,  and  torque 
coefficients,  Nusselt  number,  and  Sherwood  number  axe  determined.  Time-averaged  val¬ 
ues  of  these  fluctuating  quantities  are  also  determined.  Such  quantities  should  be  useful 
in  modelling  droplet  dispersion  and  modifications  of  heating  and  vaporization  rates  due 
to  turbulence. 
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Ill  SUMMARY  OF  RESEARCH 


1.  Unsteady  flow  interactions  between  an  advected  cylindrical 
vortex  tube  and  a  sphere 

As  a  preliminary  step  towards  understanding  the  interactions  between  a  droplet  and 
the  carrier  turbulent  flow,  the  unsteady,  three-dimensional,  incompressible,  viscous  flow 
interactions  between  a  vortical  (initially  cylindrical)  structure  advected  by  a  uniform  free 
stream  and  a  spherical  particle  held  fixed  in  space  was  investigated  numerically  for  a  range 
of  particle  Reynolds  numbers  between  20  and  100.  The  counter-clockwise  rotating  vortex 
tube  was  initially  located  ten  sphere  radii  upstream  from  the  sphere  center. 

The  computed  velocity  and  pressure  fields  provided  the  lift,  moment,  and  drag  co¬ 
efficients  on  the  sphere  as  a  function  of  time  for  a  range  of  offset  distance,  vortex  core 
radius,  and  maximum  fluctuation  velocity  induced  by  the  vortex  tube.  Initially,  the  lift 
forces  were  positive  due  to  upwash  on  the  sphere,  then  became  negative  due  to  downwash 
and  higher  fluid  velocity  near  the  bottom  of  the  sphere  when  the  vortex  tube  passed  the 
sphere. 

Varying  the  vortex  core  radius  showed  that  the  maximum  positive  and  negative  lift 
coefficients  and  the  rms  lift  coefficient  were  linearly  proportional  to  the  circulation  of  the 
vortex  tube  at  small  values  of  the  core  radius  while  they  were  linearly  proportional  only 
to  the  maximum  fluctuation  velocity  and  independent  of  the  core  radius  at  large  values 
of  the  core  radius.  For  mid-range  values  of  the  core  radius,  they  depended  on  both  the 
core  radius  and  the  maximum  fluctuation  velocity  (or  equivalently  both  the  core  radius 
and  circulation).  Some  interesting  unsteady  flow  phenomena  occurred  in  the  near  wake 
on  the  passage  of  the  vortex  tube. 

The  paper  [1]  describing  the  details  of  this  work  was  published  in  the  Journal  of  Fluid 
Mechanics  and  was  appended  to  the  September  1995  Progress  Report. 


2.  Unsteady  flow  interactions  between  a  pair  of  advected  vortex 
tubes  and  a  sphere 

The  interactions  between  a  pair  of  vortex  tubes  and  a  sphere  were  studied  in  order  to 
generalize  the  findings  from  the  previous  investigation. 

When  the  top  and  bottom  (see  Figure  2  of  Reference  2)  vortex  tubes  have  positive  and 
negative  circulations,  respectively,  the  magnitude  of  the  induced  velocity  due  to  the  vortex 
tubes  is  added  to  the  base  flow  velocity  along  the  stagnation  streamline.  This  causes  the 
pressure  at  the  stagnation  point  and  the  shear  stresses  in  the  upper  and  lower  left  regions 
to  be  higher  than  those  of  the  axisymmetric  flow  past  a  sphere,  thus  increasing  the  drag. 
On  the  other  hand,  when  the  top  and  bottom  vortex  tubes  have  negative  and  positive 
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circulations,  respectively,  the  induced  velocity  due  to  the  vortex  tubes  is  subtracted  from 
the  base  flow  velocity  along  the  stagnation  streamline.  This  causes  the  pressure  at  the 
stagnation  point  and  the  shear  stresses  in  the  upper  and  lower  left  regions  to  be  lower 
than  those  of  the  axisymmetric  flow  past  a  sphere,  thus  reducing  the  drag.  The  lift  and 
moment  are  zero  for  this  symmetric  configuration. _  .  — . 

The  interactions  between  a  sphere  and^fike-rotatin/a  pair  of  cylindrical  vortex  tubes 
initially  located  ten  radii  upstream  from  the  center  of  the  sphere  were  investigated.  The 
lift  and  moment  coefficients  of  the  sphere  interacting  with  a  pair  of  vortex  tubes  as  a 
function  of  time  are  nearly  identical,  respectively,  to  those  of  the  sphere  interacting  with 
a  single  vortex  tube  if  the  separation  distance  between  the  tube  centers  is  less  than  2 
vortex  tube  diameter  for  the  lift  coefficient  and  less  than  y/a  vortex  tube  diameter  for 
the  moment  coefficient;  here,  vmaxt  instead  of  vmax  is  used  in  the  case  of  a  pair  of  vortex 
tubes,  where  vmax  is  the  maximum  induced  velocity  due  to  one  vortex  without  presence 
of  the  other,  vmaxt  is  the  total  maximum  induced  velocity  due  to  the  pair  of  vortices  and 
a  is  the  radius  of  initial  vortex  core  normalized  by  the  sphere  radius.  In  particular,  lift 
and  moment  coefficients  are  linearly  proportional  to  the  maximum  induced  velocity.  The 
moment  coefficient  is  negligible  compared  to  the  lift  coefficient. 

The  paper  [2]  describing  this  work  in  detail  was  submitted  to  the  International  Journal 
of  Multiphase  Flow  and  is  appended  to  this  report. 


3.  Interactions  of  an  array  of  vortex  tubes  of  like  rotation  with 
a  moving  sphere 

The  two-dimensional  trajectories  of  a  spherical  particle  interacting  with  an  array  of 
vortices  whose  sizes  are  comparable  to  the  sphere  size  were  examined.  The  time-dependent 
drag  and  lift  forces  (Ref.  [1])  for  the  case  of  a  sphere  interacting  with  a  single  vortex  were 
used  to  calculate  the  two-dimensional  trajectory  of  a  moving  spherical  particle  interacting 
with  an  array  of  vortex  tubes  of  like  rotation.  The  results  show  that  the  shear  flow  across 
the  sphere  induced  by  a  vortex  tube  is  responsible  for  the  net  deflection  of  a  sphere 
interacting  with  an  array  of  vortex  tubes.  Thus,  the  sphere  eventually  deflects  in  the 
direction  of  increasing  relative  velocity.  The  deflection  ratio  (ratio  of  sphere  final  location 
in  z  and  x  directions)  of  the  sphere  increases  with  decreasing  initial  Reynolds  number 
and  with  decreasing  density  ratio.  However,  the  total  deflection  increases  with  increasing 
the  initial  Reynolds  number  and  the  density  ratio  because  higher  momentum  causes  the 
sphere  to  travel  farther.  See  Ref.  [2]  for  details. 
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4.  The  motion  of  a  sphere  in  unsteady  axisymmetric  flows  at 
moderate  Reynolds  numbers 


A  turbulent  flow  possesses  a  wide  spectrum  of  eddy  sizes.  In  order  to  enhance  the 
understanding  of  droplet  motion  in  a  turbulent  flow,  we  have  been  investigating  the  motion 
of  a  freely  moving  sphere  interacting  with  a  large  vortex  tube  whose  size  is  of  the  order 
of  an  integral  length  scale,  i.e.  at  the  other  end  of  the  spectrum  relative  to  the  previous 
cases.  The  goals  of  this  study  are  to  derive  not  only  a  relationship  of  interaction  between 
a  small  droplet  and  a  large  vortex  but  also  an  accurate  equation  for  the  motion  of  a 
spherical  droplet  in  finite-Reynolds-number  turbulent  flows. 

The  equations  of  sphere  motion  proposed  by  previous  workers  were  examined  and 
compared  with  the  results  of  the  full  Navier-Stokes  equations  for  unsteady,  axisymmetric 
flow  around  a  freely  moving  sphere  initially  injected  into  a  stationary  or  oscillating  fluid. 
As  a  result,  we  have  proposed  a  modified  equation  of  sphere  motion  and  demonstrated  its 
superiority  to  the  previously  proposed  equations  of  sphere  motion.  This  work  is  described 
in  the  AIAA  preprint  96-0081  by  Kim,  Elghobashi  &  Sirignano  [3]  which  is  an  addendum 
to  this  report. 

The  new  equation  contains  a  modified  history  term  with  an  integral  kernel  weighted  by 
the  acceleration  magnitude  and  a  new  term  accounting  for  the  initial  velocity  difference 
between  the  particle  and  the  carrier  fluid.  The  weighting  function  contains  the  time 
derivative  of  the  relative  velocity  Max  and  the  ratio  4>r  of  Mai  to  Max-  Mai,Ma2,  and  <j>r 
are  defined  by 


Max 


2  a 

|u  —  v\2 


d\u  —  v 
dt 


MA2(t )  = 


(2a)2 

|«  —  d|3 


<P\u  —  v 
dt2 


<t>T  (t)  = 


MA2 

Max 


These  dimensionless  groups  can  be  introduced  through  dimensional  analysis  to  obtain  the 
forces  on  the  particle  of  unsteady  motion. 

The  new  proposed  equation  of  particle  motion  is  expressed  as: 


dv 


Du  dv ' 
~Dt  ~  ~dt , 


=  | CDatdTra2pj  |  u  -  v  |  (u  -  v)  +  ^ 

67 vfifd  j  K(t  —  r>r)~“j — —dr  +  fa/ifa  Ki(t)  [u(0)  —  u(0)] 
with  the  integral  kernel  K(t  —  r,  r)  given  by 


+  m/~Dt  +  (mP_m/)^  + 


(la) 
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K{t  -  r,  r )  = 
G(t)  = 

P  = 


ir(t  —  t)vj 


a i 


0.5/ci 


+  G{t) 


7 r  |u(r)  —  v(r)|: 


2  avfpH(Rct) 


(t  -  rf 


1  1/cl 


-Cl 


1  +  l(r) 

C2 


1  +  + 
fH  =  0.75  +  c5Ret(r ) , 
and  with  the  function  iifi(t)  given  by 


(lb) 

(lc) 

(ld) 

(le) 


Ki  (t) 
Gi 


x  |u(0)  —  n(0)[3  21 1/ci  ) 

2  avff%(Reto)  J  J 


1  +  c^Re^'25 (pr  +  0.5)  0-5 


(lf) 

(lg) 


where  Ret  =  |n(r)  —  u(r)|2a/i//  and  Ret0  =  |u(0)  —  v(0)\2a/vj. 

The  six  constants  c,  (i  =  1,  ..,6)  in  the  above  equations  are  determined  by  comparing 
the  solutions  from  Equation  (la)  with  those  from  the  Navier-Stokes  equations  and  given 
for  0  <  a/  <  oo  by 


Ci  =  2.5,  C5  =  0.126,  C6  =  15, 

c2  =  45.5,  c3  =  0.03,  c4  =  0.1 .  (lh) 

u'  is  the  dimensionless  frequency  of  the  flow  and  defined  by  u'  =  u>a/v0 ,  where  v0  is  the 
initial  injection  velocity  of  the  particle. 

For  low  frequencies  0  <  w'  <  0.3,  the  following  values  of  c2,  c3,  and  c4  provide  slightly 
better  results  than  those  of  Equation  (lh).  The  values  of  ci,  C5,  and  c$  are  kept  fixed  as 
in  (lh). 

c2  =  13.9,  c3  =  0.12,  c4  =  0.5  .  (li) 


5.  The  motion  of  a  sphere  in  unsteady  three-dimensional  flows 
at  moderate  Reynolds  numbers 

In  order  to  extend  Equation  (la)  to  a  vector  form  for  two-  or  three-dimensional  case, 
we  consider  the  time-dependent,  three-dimensional,  incompressible,  flow  around  a  small 
spherical  solid  particle  injected  into  a  counter-clockwise  rotating  large  vortex  which  is 
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located  at  the  origin  of  the  coordinates  (X,  Y,Z)  fixed  in  space.  The  net  gravity  force 
acting  on  the  particle  is  neglected.  The  origin  of  a  nonrotating  noninertial  reference  frame 
(x,y,z)  is  chosen  at  the  center  of  the  particle. 

All  the  variables  are  nondimensionalized  using  the  sphere  radius  a  as  the  characteristic 
length  and  the  initial  injection  velocity  of  the  sphere  vQ  as  the  characteristic  velocity. 

The  initial  velocity  field  induced  by  the  vortex  tube  is  analytically  computed  by  con¬ 
sidering  the  evolution  of  a  point  vortex  and  is  given  by 

u*  =  2^  -  exp(-i 7JTTT7))]  <2a) 

, _ ?L_ 

°  4i//(1.12)a’ 

where  R  =  \/X2  +  Y2,  <j>  =  arctany/X,  and  t0  is  a  parameter  defined  by  initial  size  of 
the  vortex  core  ( <r0 )  and  the  fluid  kinematic  viscosity  (t//). 

The  initial  maximum  induced  velocity  umaxo  due  to  the  vortex  tube  occurs  on  the 
edge  of  the  vortex  core  R  =  cr0  at  t  =  0  and  can  be  obtained  from  Equation  (2a).  The 
pressure  field  due  to  the  vortex  tube  is  also  analytically  computed  by  integrating  the 
radial  component  of  the  momentum  equation  which  is  pU^/R  =  dp/dR. 

It  is  assumed  that  the  vortex  is  so  large  that  the  flow  field  induced  by  the  vortex  is 
not  affected  by  the  moving  sphere  except  locally  in  the  region  around  the  sphere.  The 
boundary  conditions  of  the  flow  field  are  obtained  by  superimposing  the  sphere  velocity 
and  the  induced  velocity  due  to  the  vortex  tube  at  the  computational  outer  boundary 
relative  to  the  sphere,  and  the  pressure  field  due  to  the  vortex  tube  is  also  imposed  at  the 
computational  outer  boundary.  Although  the  flow  computation  is  three-dimensional,  the 
particle  path  remains  in  the  plane  of  symmetry  (X  —  Y  plane).  Therefore,  the  trajectory 
computation  is  made  for  the  case  of  a  freely  moving  particle  in  two-dimensions.  After 
computing  the  forces  on  the  sphere,  the  deceleration  (or  acceleration)  of  the  sphere  is 
obtained  via  Newton’s  second  law  of  motion,  and  then  the  new  location  of  the  sphere  is 
obtained. 

Preliminary  results  of  the  large  vortex  study  show  that  the  lift  force  on  the  particle 
is  much  smaller  than  that  due  to  the  small  vortices  with  the  same  maximum  induced 
velocity.  Our  direct  solution  of  the  three-dimensional  Navier-Stokes  equations  over  a 
freely  moving  particle  confirms  that  due  mainly  to  the  drag  force,  the  particle  travels 
in  a  curved  trajectory  that  depends  on  the  direction  of  vortex  rotation  and  the  Stokes 
number  of  the  particle.  This  indicates  that  the  mechanism  of  particle  dispersion  due  to 
the  interaction  with  small  vortices  is  quite  different  from  that  due  to  the  interaction  with 
a  large  vortex.  Refer  to  [3]  for  details  of  the  work. 


(2b) 
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6.  Unsteady  thermal  interactions  between  a  liquid  sphere  and 
an  advecting  vortex 

Within  the  scope  of  our  planned  study  on  the  thermal  aspects  of  the  droplet-turbulence 
interaction,  we  have  completed  quantifying  the  effect  of  an  advecting  vortex  near  the 
droplet  on  the  droplet  convective  heat  transfer.  Also,  prompted  by  recent  publications, 
we  have  included  the  effect  of  Marangoni  stresses  in  the  dynamics.  We  intend  to  study 
in  the  future  a  coupled  effect  of  the  gas  field  temperature  stratification  and  advecting 
vortical  structures  on  the  droplet  heating. 

We  successfully  completed  our  study  of  droplet  heating  influenced  by  an  advecting 
vortex.  Since  the  problem  is  unsteady  by  nature,  we  paid  particular  attention  to  time- 
averaged  and  root-mean-squared  values  of  the  droplet  Nusselt  number.  All  the  existing 
available  correlation/jfor  droplet  heat  transfer  are  for  a  steady  axisymmetric  flow;  thus, 
they  can  not  be  applied  when  there  exists  asymmetry  or  unsteadiness  in  the  gas  field, 
such  as  that  induced  by  an  advecting  vortex.  It  is  noteworthy  that,  in  practice,  droplets 
do  not  experience  a  perfectly  symmetric  field  but  an  asymmetric  one.  This  introduces 
serious  limitation  on  the  applicability  of  existing  correlations. 

We  have  successfully  produced  the  first  known  self-similar  correlation  predicting  the 
unsteady  droplet  heating  due  to  the  advecting  vortex;  this  finding  well  compliments  the 
existing  correlations  for  droplet  heating  in  symmetric  flows.  Our  results  include  several 
interesting  observations:  first,  self- similarity  in  the  droplet  heating  exists  for  the  unsteady 
droplet- vortex  interaction;  second,  it  is  not  the  vortex  initial  distance  from  the  droplet 
but  its  ratio  to  the  vortex  core  size  that  plays  a  major  role  in  droplet  heating;  this  ratio 
determines  whether  the  droplet  would  be  embedded  inside  the  vortex  inner  core  in  the 
course  of  the  interaction.  Third,  if  this  ratio  is  larger  than  a  certain  value,  the  droplet 
heating  merely  depends  on  the  vortex  circulation  and,  surprisingly,  has  little  dependence 
on  the  vortex  distance  from  the  droplet.  The  role  of  the  flow  Reynolds  number  itself 
appears  to  be  crucial  as  well:  droplet  heating  shows  stronger  response  to  the  same  vortex 
in  a  flow  with  a  higher  Reynolds  number.  In  general,  it  was  found  that  the  droplet  Nusselt 
number  variation  is  qualitatively  similar  to  change  in  droplet  fluid  dynamic  properties 
observed  in  previous  studies  of  this  research  group;  coupling  between  the  fluid  dynamics 
and  heat  transfer  phenomena  justifies  this  similarity.  These  findings  [4]  were  presented  at 
the  24th  Aerospace  Sciences  Meeting,  January  1996,  in  Reno,  Nevada,  and  are  attached 
to  this  report  as  an  Addendum.  An  expanded  version  of  this  paper  is  being  submitted  to 
the  International  Journal  of  Heat  and  Mass  transfer  for  publication. 

Recent  publications  by  Niazmand  et  al  [5]  and  Shih  and  Megaridis  [6]  signify  some  im¬ 
portance  of  thermocapillary  effects  on  the  droplet  response  in  a  hot  gaseous  axisymmetric 
environment.  Most  previous  publications  in  the  field  of  spray  droplet  computations  have 
neglected  these  effects.  Following  these  findings,  we  included  the  effect  of  surface  tension, 
the  so-called  Marangoni  stresses,  in  the  force-balance  at  the  droplet-gas  interface  and  so 
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in  the  numerical  simulation.  This  task  is  completed  and  the  code  has  been  benchmarked 
after  this  upgrade.  We  have  successfully  reproduced  the  effect  of  thermocapillary  stresses 
on  the  droplet  internal  circulation  underlined  in  said  publications,  the  most  important 
ones  being  the  slow-down  on  the  droplet  internal  circulation  and  the  detachment  and 
gradual  disappearance  of  the  droplet  wake. 

Droplets  in  a  combustion  chamber  do  not  experience  a  perfectly-uniform  temperature 
field  but  a  spatially  varying  one.  The  effect  of  temperature  stratification,  coupled  with 
the  influence  of  vortical  structures,  therefore  pose  interesting  questions  on  the  dynamics 
of  the  said  coupling;  e.g.  it  is  unclear  how  the  vortex-induced  variations  in  the  droplet 
heating  [4]  will  change  due  to  gas-field  temperature  stratification.  In  spite  of  the  serious 
practiced  applications  for  such  findings,  there  exist  no  available  information  shedding  light 
on  the  coupled  temperature  stratification-vortical  stratification  effects  on  droplet  heating. 

We  therefore  intend  to  study  simultaneous  thermal  and  aerodynamic  effects  of  the 
temperature  field  and  the  vortical  structures  on  the  droplet  heating  in  a  combustion 
chamber.  We  assume  a  profile  for  the  surrounding  gas  temperature  whose  length-scale  of 
its  gradient  is  comparable  to  that  of  the  droplet  and  the  vortex.  Otherwise  stated,  there 
are  three  important  length  scales  in  this  study:  the  vortex  core  size  ,  droplet  size,  and 
the  length  scale  for  the  temperature  stratification.  Special  interest  occurs  when  the  three 
length  scales  are  comparable. 
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ADDENDUM  TO  SECTIONS  2  and  3 


UNSTEADY  FLOW  INTERACTIONS  BETWEEN 
A  PAIR  OF  ADVECTED  VORTEX  TUBES  AND 
A  RIGID  SPHERE 

by 

Inchul  Kim,  Said  Elghobashi,  and  William  A.  Sirignano 

Department  of  Mechanical  and  Aerospace  Engineering 
University  of  California,  Irvine 
Irvine,  CA  92717 


Abstract 

An  idealized  representation  of  the  interaction  of  spherical  particles  with  turbulent  eddies  of 
comparable  length  scale  is  considered  by  means  of  a  three-dimensional,  unsteady  finite-difference 
Navier-Stokes  solution  of  the  interaction  between  a  fixed  rigid  sphere  and  a  pair  ol  advecting 
vortex  tubes.  First,  a  doubly  symmetric  interaction  with  vortices  of  opposite  rotation  is  consid¬ 
ered.  The  resulting  time-dependent  drag  differs  from  the  drag  in  axisymmetric  flows;  however, 
the  lift  and  torque  on  the  sphere  remain  zero.  Next,  an  interaction  with  two  vortices  of  like 
rotation  is  studied.  Here,  non-zero  lift  and  torque,  as  well  as  drag  deviation  from  the  axisym¬ 
metric  case  occur  and  would  result  in  a  deflection  in  the  trajectory  of  a  nonfixed  sphere.  The 
flow  in  this  case  behaves  like  that  of  a  single  vortex.  Finally,  a  linear  array  of  like-rotating 
vortices,  interacting  with  a  freely  moving  sphere,  is  considered.  The  two-dimensional  deflection 
depends  strongly  upon  the  sphere /fluid  density  ratio  and  initial  sphere  Reynolds  number.  Lift 
and  moment  coefficients  are  shown  to  be  linearly  proportional  to  the  maximum  induced  velocity 
due  to  the  vortices.  Moment  coefficients  are  an  order  of  magnitude  less  than  lift  coefficients. 

Key  Words:  unsteady  flow  over  a  sphere,  sphere- vortex  interaction 
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1  Introduction 


The  interactions  between  vortical  structures  and  spherical  particles  or  droplets  is  of 
primary  practical  interest  in  many  particle-laden  flows.  These  interactions  modify  the 
trajectories  of  individual  particles  and  cause  dispersion  in  a  spray  or  cloud.  Also,  they 
can  modify  heat  and  mass  transfer  rates  for  the  particles.  There  has  been  long-term 
interest  in  the  effects  of  turbulent  eddies  which  contain  most  the  energy  and  whose  sizes 
are  orders  of  magnitude  larger  than  the  particle  diameters. 

However,  a  need  exists  to  study  the  interactions  of  a  particle  with  vortical  structures 
that  are  smaller  than,  comparable  in  size  to,  or  only  a  few  times  larger  than  the  sphere. 
These  small  structures  have  the  potential  to  produce  the  largest  modifications  to  the 
boundary  layer  and  near  wake  of  the  sphere. 

In  this  paper,  an  idealized  representation  of  those  interactions  is  made  by  considering 
the  viscous,  incompressible,  unsteady,  three-dimensional  flow  associated  with  a  pair  of 
initially  cylindrical  vortex  tubes  advecting  past  the  sphere.  This  study  builds  upon  the 
previous  study  of  Kim,  Elghobashi  k  Sirignano  (1995),  hereafter  identified  as  KES.  In 
that  paper,  the  authors  examined  the  unsteady,  three-dimensional  interactions  between 
a  single  advected  cylindrical  vortex  tube  and  a  fixed  spherical  particle  whose  diameter  is 
of  the  same  order  of  magnitude  as  the  initial  diameter  of  the  vortex.  That  study  served 
as  a  first  step  towards  better  understanding  the  two-way  interactions  between  small-scale 
turbulence  and  the  particle.  Here  we  extend  this  work,  with  the  same  goal,  to  study  the 
interactions  between  a  pair  of  advected  vortex  tubes  and  a  stationary  spherical  particle. 

In  the  earlier  study  (KES)  the  particle  Reynolds  number  based  on  the  freestream 
velocity  and  the  particle  diameter  was  in  the  range  20  <  Re  <  100.  The  initial  size  of 
the  cylindrical  vortex  tube  was  in  the  range  0.25  <  a  <  4,  where  a  is  the  radius  of  the 
vortex  tube  normalized  by  that  of  the  particle.  We  found  that  the  maximum  positive 
lift  coefficient  and  the  rms  lift  coefficient  of  the  sphere  are  linearly  proportional  to  the 
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circulation  of  the  vortex  tube  at  small  values  of  a.  However,  at  large  values  of  <7,  they 
are  linearly  proportional  to  the  maximum  induced  velocity  due  to  the  vortex  tube  but 
independent  of  a. 

In  the  present  paper,  both  like-rotating  and  opposite-rotating  vortical  pairs  axe  con¬ 
sidered.  In  the  opposite- rotating  case,  only  a  symmetric  configuration  (see  figure  1)  is 
examined.  Asymmetric  configurations  with  opposite  rotations  are  left  for  future  studies. 
More  attention  is  given  to  like-rotating  pairs  because  they  have  the  greatest  effect  on  lift 
and  torque.  We  expect,  therefore,  that  deflections  in  the  trajectories  and  dispersion  of 
sprays  and  clouds  will  be  greater  in  this  case  of  like  rotation.  The  case  of  a  "train”  of 
vortices  advecting  past  the  sphere  at  prescribed  intervals  is  also  examined. 

Our  specific  objectives  are  to  study: 

1.  the  detailed  flow  field  behavior  during  interaction  of  a  pair  of  vortex  tubes  with  each 
other  and  with  the  sphere, 

2.  the  relationship  between  the  lift  coefficient  of  the  sphere  and  the  maximum  induced 
velocity  due  to  the  two  vortex  tubes, 

3.  the  modification  of  the  drag  force  caused  by  the  interactions, 

4.  the  effects  of  Reynolds  number,  vortex  size,  and  initial  offset  distance  of  the  vortex, 
and 

5.  the  sphere  deflection  caused  by  the  interaction  with  the  vortices. 

The  detailed  study  of  the  interactions  between  the  particle  and  the  unsteady  veloc¬ 
ity  field  provides  fundamental  information  about  the  flow  behavior  that  can  be  used  in 
developing  mathematical  models  for  particle-laden  flows.  The  next  section  provides  a 
mathematical  description  of  the  flow  considered,  the  governing  equations,  and  the  numer¬ 
ical  solution  procedure.  Section  3  discusses  the  results  including  the  numerical  accuracy 
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issues,  the  effects  of  varying  the  parameters  listed  above,  and  the  trajectories  of  a  moving 
spherical  particle  interacting  with  an  array  of  vortex  tubes  of  like  rotation  as  an  extension 
of  the  results  of  KES  and  the  present  study.  Section  4  provides  a  summary  and  concluding 
remarks. 


2  Problem  statement  and  formulation 

2.1  Flow  description 

Consider  the  time-dependent,  three-dimensional,  incompressible,  viscous  flow  interac¬ 
tions  between  a  pair  of  symmetric,  initially  cylindrical  vortex  tubes  and  a  solid  sphere. 
The  vortex  tubes  are  moving  with  the  laminar  free  stream,  and  a  sphere  is  suddenly  placed 
and  held  fixed  in  space  as  shown  in  figure  1.  The  initial  offset  distance,  <//>  denotes  the 
shortest  distance  between  the  initial  vortical  axis  and  the  y-z  plane  which  is  parallel  to 
the  free  stream.  All  the  variables  are  nondimensionalized  using  the  sphere  radius  a'  as 
the  characteristic  length  and  as  the  characteristic  velocity,  where  the  superscript  1 
denotes  dimensional  quantity.  The  two  vortex  tubes,  having  equal  diameters  of  the  order 
of  the  sphere  diameter,  are  initially  located  ten  sphere-radii  upstream  from  the  center 
of  the  sphere.  The  effects  of  the  vortex  tubes  on  the  sphere  are  negligible  at  this  initial 
distance  because  the  magnitude  of  the  initial  velocity  field  induced  by  the  vortex  tubes 
is  less  than  2  percent  of  the  free  stream  velocity.  Far  upstream,  the  flow  is  uniform  with 
constant  velocity  U'^k  parallel  to  the  y-z  plane.  There  is  one  symmetry  plane,  the  x-z 
plane,  as  seen  in  figure  1.  A  second  symmetry  plane  (y-z)  exists  only  when  the  two  vor¬ 
tices  have  opposite  rotations.  Our  general  formulation  does  not  take  advantage  of  this 
second  symmetry. 

Note  that,  later  in  section  3.4,  the  fixed  sphere  results  will  be  employed  in  a  moving 
sphere  trajectory  analysis. 


4 


Two  coordinate  systems  are  used  in  our  formulation  following  KES:  the  Cartesian 
coordinates  (x,y,z)  and  the  nonorthogonal  generalized  coordinates  (£,  77,  ().  The  origin 
of  the  former  coincides  with  the  sphere  center.  £  is  the  radial,  77  is  the  angular,  and  ( 
is  the  azimuthal  coordinates.  The  nonorthogonal  generalized  coordinate  system  can  be 
easily  adapted  to  three-dimensional  arbitrary  geometries.  In  the  present  study,  a  spherical 
domain  is  used,  and  the  grid  reduces  to  an  orthogonal,  spherical  grid.  The  grids  are  denser 
near  the  surface  of  the  spherical  particle,  and  the  grid  density  in  the  radial  direction  is 
controlled  by  the  stretching  function  developed  by  Vinokur  (1983).  Due  to  symmetry,  the 
physical  domain  is  reduced  to  a  half  spherical  space.  The  domain  of  the  flow  is  bounded 
by  1  <  £  <  N1}  1  <  r/  <  N2,  1  <  (  <  N3,  where  £  =  1  and  N\  correspond,  respectively,  to 
the  sphere  surface  and  the  farfield  boundary  surrounding  the  sphere;  77  =  1  and  7V2  denote, 
respectively,  the  positive  z-axis  (downstream)  and  the  negative  z-axis  (upstream);  (  =  l 
and  Nz  refer,  respectively,  to  the  x-z  plane  in  the  positive  x-direction  and  the  x-z  plane 
in  the  negative  x-direction.  Uniform  spacing  (6£  =  6rj  =  6(  =  1)  is  used,  for  convenience, 
for  the  generalized  coordinates. 

The  initial  vortex  tubes  have  a  small  core  region  with  a  radius  a  (normalized  by  the 

sphere  radius).  This  core  is  defined  such  that  the  initial  velocity  induced  by  the  vortex 

tube  approaches  zero  as  the  distance  from  the  center  of  the  vortex  tube  goes  to  zero,  and 

at  distances  much  greater  than  cr,  the  induced  velocity  approaches  that  of  a  point  vortex. 

We  use  the  vortex  tube  construction  of  Spalart  (1982),  which  has  the  following  stream 

function:  r 

7/>„(x,  z,t  =  0)  =  ~2^/n[(x  -  Xj  f  +  (z  -  Zj )2  +  a2]  , 

where  T,-  is  the  nondimensional  circulation  around  the  vortex  tube  at  radius  a  and  at  the 
initial  time.  Tj  is  positive  when  the  vortex  tube  rotates  counterclockwise,  and  x}  and  z3 
denote  the  location  of  the  center  of  the  vortex  tube.  The  circulation  around  a  circular 
path  far  away  from  the  center  of  the  vortex  is  given  by  D,*  =  2Y 3.  Each  vortex  tube  can 
be  viewed  as  an  evolution  from  the  point  vortex  due  to  the  cylindrical  viscous  diffusion. 
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(1) 


The  stream  function  for  a  pair  of  vortex  tubes  is  given  by 

ij)v{x,z,i  =  0)  =  -£  -  Xj)2  +  (z  -  Zjf  +  <r2] 

;=i  27r 

2.2  Governing  equations  and  boundary  conditions 

The  continuity  and  momentum  equations  and  the  initial  and  boundary  conditions  are 
nondimensionalized  using  the  sphere  radius  a'0  as  the  characteristic  length  and  U1^  as  the 
characteristic  velocity. 


W  =  0  (2a) 

^  +  v-yy  =  -Vp  +  -|-v2y  (2b) 

at  tie 

The  governing  equations  (2a)  and  (2b)  are  cast  in  terms  of  the  generalized  coordinates 
(^^0  to  treat  a  three-dimensional  body  of  arbitrary  shape.  The  numerical  integration 
is  performed  using  a  cubic  computational  mesh  with  equal  spacing  (6£  =  Srj  =  =  1). 

The  velocities  on  the  sphere  surface  are  zero  due  to  the  no-slip  condition,  and  the 
pressure  on  the  sphere  is  obtained  from  the  momentum  equation.  The  detailed  equations 
describing  the  boundary  and  initial  conditions  are  given  in  KES  and  thus  will  not  be 
repeated  here.  The  only  difference  from  KES  is  that  the  initial  pressure  is  estimated  as 
zero  over  the  whole  computational  domain.  This  estimation  is  corrected  by  the  pressure 
correction  equation  and  iteration  procedure  (see  section  2.3  for  details). 

The  only  nondimensional  groupings  appearing  in  the  equations  and  initial  and  bound¬ 
ary  constraints  are  the  sphere  Reynolds  number,  vortex  tube  radius,  offset  distance,  and 
vortex  circulation  (or  vortex  Reynolds  number). 

The  equations  evaluating  the  drag,  lift,  and  moment  coefficients  are  given  in  KES 
and  thus  will  not  be  repeated  here.  The  lift  force  is  assumed  positive  when  it  is  directed 
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toward  the  positive  x-axis.  Due  to  symmetry,  only  the  y-component  of  the  moment  is 
non-zero  and  is  assumed  positive  in  the  counter-clockwise  direction. 


2.3  Numerical  solution 

A  three-dimensional,  implicit,  finite-difference  algorithm  has  been  developed  to  solve 
simultaneously  the  set  of  the  discretized  partial  differential  equations.  The  method  is 
based  on  an  Altemating-Direction-Predictor-Corrector  (ADPC)  scheme  to  solve  the  time- 
dependent  Navier-Stokes  equations.  ADPC  is  a  slight  variation  of  Altemating-Direction- 
Implicit  (ADI)  method  and  implemented  easily  when  embedded  in  a  large  iteration  scheme 
(Patnaik  1986,  Patnaik  et  al.  1986).  The  control  volume  formulation  is  used  to  develop  the 
finite-difference  equations  from  the  governing  equations  with  respect  to  the  generalized 
coordinates  (£,r?,C).  One  of  the  advantages  of  the  control  volume  formulation  is  that 
mass  and  momentum  are  conserved  over  a  single  control  volume,  and  hence  the  whole 
domain  regardless  of  the  grid  fineness.  An  important  part  of  solving  the  Navier-Stokes 
equations  in  primitive  variables  is  the  calculation  of  the  pressure  field.  In  the  present  work, 
a  pressure  correction  equation  is  employed  to  satisfy  indirectly  the  continuity  equation 
(Anderson  et  al.  1984).  The  pressure  correction  equation  is  of  the  Poisson  type  and  is 
solved  by  the  Successive-Over-Relaxation  (SOR)  method. 

The  overall  solution  procedure  is  based  on  a  cyclic  series  of  guess-and-correct  opera¬ 
tions.  The  velocity  components  are  first  calculated  from  the  momentum  equations  using 
the  ADPC  method,  where  the  pressure  field  at  the  previous  time  step  is  employed.  This 
estimate  improves  as  the  overall  iteration  continues.  The  pressure  correction  is  calculated 
from  the  pressure  correction  equation  using  the  SOR  method,  and  new  estimates  for  pres¬ 
sure  and  velocities  are  obtained.  This  process  continues  until  the  solution  converges  at 
each  time  step. 
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3  Results  and  discussion 


In  subsections  (3.1)  and  (3.2),  we  discuss  the  three-dimensional  interactions  of  a 
sphere  and  a  pair  of  vortex  tubes  of  opposite  rotation.  In  subsection  (3.3),  we  examine 
the  three-dimensional  interactions  of  a  sphere  and  a  pair  of  vortex  tubes  of  like  rotation. 
In  subsection  (3.4),  we  investigate  the  trajectories  of  a  moving  sphere  interacting  with  an 
array  of  vortex  tubes  of  like  rotation. 

Testing  the  accuracy  of  our  numerical  solution  has  been  performed  and  discussed 
earlier  in  KES.  The  51  x  51  x  51  grid  is  used  in  the  following  calculations.  The  run  for 
the  interaction  between  a  single  vortex  tube  and  a  sphere  at  Reynolds  number  100  with 
the  51  x  51  x  51  grid  required  4.95  mega  words,  a  dimensionless  time  step  of  At  =  0.002, 
and  a  total  time  of  4  cpu  hours  on  Cray  C-90  for  the  final  time  of  tj  =  24.5.  Each  time 
step  takes  about  1.18  cpu  seconds. 

3.1  Interactions  of  a  sphere  and  a  pair  of  vortex  tubes  with 
top-positive  and  bottom-negative  circulations 

We  consider  the  interactions  of  a  pair  of  vortex  tubes  advected  by  the  free  stream  and 
a  sphere  suddenly  placed  in  the  flow  and  held  fixed  in  space.  The  two  cylindrical  vortex 
tubes  are  initially  of  the  same  size  and  rotating  opposite  to  each  other  with  top-positive 
and  bottom-negative  circulations  as  shown  in  figure  1.  The  y-z  plane  is  half  way  between 
the  two  tubes  so  that  the  offset  distance  of  one  vortex  tube  is  the  negative  of  the  offset 
distance  of  the  other.  The  center  of  the  each  vortex  tube  is  located  at  10  sphere-radii 
upstream  from  x-y  plane  containing  the  center  of  the  sphere.  The  base  case  calculation 
is  that  of  Re  —  100,  d0jf  =  ±1.5,  and  a  =  1. 

Initially,  each  vortex  tube  has  its  maximum  induced  velocity  vmax  located  at  the  edge 
of  the  core.  Because  the  velocity  and  vorticity  fields  induced  by  one  vortex  tube  influence 
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those  by  the  other,  the  total  maximum  induced  velocity,  vmaxt,  due  to  the  two  vortex  tubes 
depends  on  their  size  and  separation  distance  and  is  in  the  range  0  <  vmaxt  <  2vmax.  The 
total  maximum  induced  velocity  vmaxt  equals  zero  when  \d0fj\  =  0,  2  vmax  when  \d0jf\  =  a, 
and  vmax  when  \d0jj\  »  1.  For  example,  vmaxi  is  0.738  for  vmax  =  0.4,  doff  =  ±1.5,  and 
a  =  1.  The  base  case  calculation  is  that  of  Re  =  100,  d0fj  =  ±1.5,  and  <7  =  1.  Note  that 
the  lift  and  torque  on  the  sphere  are  zero  due  to  the  flow  symmetry  in  upper  and  lower 
regions  of  the  sphere 

In  order  to  describe  the  flow  structure,  we  first  consider  the  pseudo-streamlines  and 
vorticity  contours  in  the  x-z  symmetry  plane,  defined  as  the  principal  plane ,  where  the 
strongest  interactions  occur  between  the  vortical  structures  and  the  sphere.  The  line 
connecting  the  front  and  rear  stagnation  points  in  the  standard  axisymmetric  flow  over 
a  single  sphere,  which  is  the  x  =  0  line  in  the  principal  plane,  will  be  used  a s  a  reference 
line.  We  refer  to  the  region  above  that  line  as  the  ‘upper’  region  and  that  below  the  line 
as  the  ‘lower’  region. 

The  pseudo-streamlines  are  obtained  from  the  pseudo-stream  function  which  is  defined 
by  assuming  that  the  velocity  field  in  the  principal  plane  does  not  change  in  the  direction 
normal  to  that  plane  and  by  using  the  two-dimensional  stream  function  definition.  The 
sphere  surface  in  the  principal  plane  is  used  as  a  reference  streamline  (ipps  =  0).  We  note 
that  a  real  stream  function  ip  cannot  be  defined  and  calculated  from  the  velocity  in  the 
principal  plane  due  to  the  existence  of  a  divergence  associated  with  the  third  component 
of  velocity.  Nevertheless,  for  descriptive  purposes  only,  it  is  convenient  to  use  the  two- 
dimensional  stream  function  definition  to  present  descriptions  of  the  flow  pattern. 

Figures  2(a)- (f)  display  the  pseudo-streamlines  (left  column)  and  the  contour  lines  of 
y-component  vorticity  (right  column)  in  the  principal  plane  at  t  =  0,  3,  6,  9,  12,  and  15 
for  Re  =  100,  d0jj  =  ±1.5,  a  =  1  with  vmaxt  =  0.738  (vmax  =  0.4).  The  contour  values  of 
the  pseudo-streamlines  are  0,  ±0.02,  ±0.1,  ±0.3.  The  contour  values  of  the  vorticity  are 
±0.4,  ±0.8,  ±1.4,  ±2,  with  the  highest  magnitude  at  the  sphere  surface. 
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Comparing  the  vorticity  contours  in  figures  2(a)-(f)  with  those  of  the  single  vortex 
tube  in  figure  5  of  KES,  we  see  that  the  two  vortex  tubes  move  downstream  faster  than 
the  single  vortex.  This  additional  acceleration  occurs  because  the  velocity  magnitude  at 
the  center  of  each  vortex  tube  equals  that  of  the  base  flow  plus  that  induced  by  the  other 
vortex  tube. 

The  distance  between  the  top  pseudo-streamline  and  the  bottom  pseudo-streamline  in 
figures  2(a)-(c)  is  narrower  on  the  segment  connecting  the  vortex  tube  centers  than  any 
other  place  along  the  stagnation  pseudo-streamline.  This  indicates  that  the  velocity  near 
the  middle  of  the  segment  between  the  vortex  tube  centers  is  higher  than  any  other  place 
along  the  stagnation  pseudo-streamline.  The  induced  velocity  due  to  the  vortex  tubes  is 
added  to  the  base  flow  near  the  stagnation  pseudo-streamline. 

Figure  3  shows  the  drag  coefficients  of  the  sphere  as  a  function  of  time  for  the  same 
parameters  as  above.  The  drag  coefficients  are  obtained  with  four  different  total  maximum 
induced  velocities  due  to  the  vortex  tubes,  vmaxt  =  0.185,  0.369,  0.554,  and  0.738  (vmax  = 
0.1,  0.2,  0.3,  and  0.4).  The  temporal  behavior  of  the  drag  coefficients  is  different  from  that 
of  the  case  of  the  pair  of  vortex  tubes  of  like  rotation  as  will  be  shown  in  section  3.3.  The 
time-averaged  value  of  the  deviation  of  the  drag  coefficient  from  that  of  the  axisymmetric 
flow  past  a  sphere  for  all  values  of  vmaxi  is  not  negligible  and  increases  linearly  with  vmaxi. 
The  time-averaged  drag  coefficient  Cn,ave  may  be  expressed  by 

=  Cu.oxt  "l"  Vmaxt  5  (3) 

where  the  constant  /?  =  0.27,  and  Co,axi  is  the  time  averaged  value  of  the  drag  coefficient 
in  the  case  of  axisymmetric  flow  {vmaxt  =  0).  The  drag  coefficients  reach  their  maximum 
at  about  t  =  9  (see  figure  3).  The  maximum  drag  coefficient  Cc)max  can  be  expressed 
approximately  by  equation  (3)  but  with  j3  =  1.05,  and  C D,ari  here  is  the  local  value  of  the 
axisymmetric  drag  coefficient  at  the  time  of  Cnimox.  Because  the  top  and  bottom  vortex 
tubes  have  positive  and  negative  circulations,  respectively,  the  induced  velocity  due  to  the 
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vortex  tubes  adds  its  magnitude  to  the  base  flow  along  the  stagnation  pseudo-streamline. 
This  increased  velocity  causes  the  pressure  at  the  stagnation  point  and  the  shear  stresses 
in  the  upper  and  lower  left  regions  to  be  higher  than  those  of  the  axisymmetric  flow  past 
a  sphere.  As  a  consequence,  the  drag  is  increased. 

3.2  Interactions  of  a  sphere  and  a  pair  of  vortex  tubes  with 
top-negative  and  bottom-positive  circulations 

We  consider  the  same  initial  flow  geometry  and  parameters  as  those  in  section  3.1  but 
for  a  pair  of  vortex  tubes  with  top-negative  and  bottom-positive  circulations. 

Figures  4(a)-(f)  display  the  pseudo-streamlines  (left  column)  and  the  contour  lines  of 
y-component  vorticity  (right  column)  in  the  principal  plane  at  t  =  0,  3,  6,  9,  12,  and  15 
for  Re  =  100,  d0jj  =  ±1.5,  <7  =  1  with  vmaxi  =  0.738  (vmax  =  0.4).  The  contour  values  of 
the  pseudo- streamlines  and  the  vorticity  are  the  same  as  those  in  the  previous  section. 

The  vorticity  contours  in  figures  4(a)-(f)  show  that  the  two  vortex  tubes  move  down¬ 
stream  slower  than  the  single  vortex  tube  in  figure  5  of  KES.  This  relative  deceleration 
occurs  because  the  velocity  magnitude  at  the  center  of  each  vortex  tube  equals  that  of 
the  base  flow  minus  that  induced  by  the  other  vortex  tube. 

The  distance  between  the  top  pseudo-streamline  and  the  bottom  pseudo-streamline 
in  figures  4(a)-(c)  is  broader  near  the  segment  connecting  the  vortex  tube  centers  than 
any  other  place  along  the  stagnation  streamline.  This  indicates  that  the  velocity  near  the 
middle  of  the  segment  between  the  vortex  tube  centers  is  lower  than  those  upstream  or 
downstream  of  the  vortex  tubes  along  the  stagnation  streamline.  The  induced  velocity 
due  to  the  vortex  tubes  is  subtracted  from  the  base  flow  near  the  stagnation  pseudo- 
streamline. 

Figure  5  shows  the  drag  coefficients  of  the  sphere  as  a  function  of  time  for  the  same 
parameters  as  used  in  section  3.3.1.  The  drag  coefficients  are  obtained  for  four  different 
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total  maximum  induced  velocities  due  to  the  vortex  tubes,  vmaxt  =  0.185,  0.369,  0.554, 
and  0.738  {vmax  =  0.1,  0.2,  0.3,  and  0.4).  The  temporal  behavior  of  the  drag  coefficients 
is  different  from  that  of  the  case  of  the  pair  of  vortex  tubes  of  like  rotation  as  will  be 
shown  in  section  3.3.  The  time-averaged  value  of  the  deviation  of  the  drag  coefficient 
from  that  of  the  axisymmetric  flow  past  a  sphere  for  all  values  of  vmaxt  is  n°t  negligible 
and  decreases  linearly  as  vmaxt  increases.  The  time-averaged  drag  coefficient  Co,ave  may 
be  expressed  by  by  equation  (3)  but  with  the  proportionality  constant  /?  =  —0.28,  and 
Cn,axi  is  the  time-averaged  value  of  the  drag  coefficient  in  the  case  of  axisymmetric  flow 
(vmaxt  =  0).  The  drag  coefficients  reach  their  minimum  at  about  t  =  11  (see  figure  5). 
The  minimum  drag  coefficient  Cn,min  can  be  expressed  approximately  by  equation  (3)  but 
with  /3  =  -0.95,  and  Cn,axi  here  is  the  local  value  of  the  axisymmetric  drag  coefficient  at 
the  time  of  Co,min-  Because  the  top  and  bottom  vortex  tubes  have  negative  and  positive 
circulations,  respectively,  the  induced  velocity  due  to  the  vortex  tubes  is  subtracted  from 
the  base  flow  velocity  along  the  stagnation  streamline.  This  causes  the  pressure  at  the 
stagnation  point  and  the  shear  stresses  in  the  upper  and  lower  left  regions  to  be  lower 
than  those  of  the  axisymmetric  flow  past  a  sphere.  Thus,  the  drag  is  reduced. 

3.3  Interactions  of  a  pair  of  vortex  tubes  of  like  rotation  and  a 
sphere 

We  consider  the  same  initial  flow  geometry  and  parameters  as  those  in  section  3.1 
but  for  a  pair  of  vortex  tubes  of  like  rotation.  The  base  case  calculation  is  that  of 
Re  =  100,  do}}  —  ±1.5,  and  <j  =  1. 

Initially  each  vortex  tube  has  its  maximum  induced  velocity  vmax  located  at  the  edge 
of  the  core.  Because  the  velocity  and  vorticity  fields  induced  by  one  vortex  tube  influence 
those  by  the  other,  the  total  maximum  induced  velocity,  umaxt,  due  to  the  two  vortex  tubes 
depends  on  their  size  and  separation  distance  and  is  in  the  range  vmax  <  vmaxt  <  2umoI. 


12 


Vmaxt  equals  2vmax  when  |<f0//|  =  0  and  equals  vmax  when  \dojj\  »  1.  For  example, 
t) maxt  “*  0*59  for  VfjmX  —  0.4,  dojf  —  ±1.5,  and  o  —  1* 

In  subsections  3.3.1  and  3.3.2,  we  investigate  the  base  case.  In  subsections  3.3.3  and 
3.3.4,  we  discuss  the  effects  of  the  size  and  the  offset  distance  of  the  vortex  tubes  and 
Reynolds  number,  respectively. 

3.3.1  Flow  structure 

Figures  6(a)-(f)  display  the  pseudo-streamlines  (left  column)  and  the  contour  lines  of 
y-component  vorticity  (right  column)  in  the  principal  plane  at  t  =  1,  6, 10,  15,  21,  and  30 
for  Re  =  100,  d0jj  =  ±1.5,  a  =  1,  and  vmaxt  =  0.59  (vmax  =  0.4).  The  contour  values  of 
the  pseudo-streamlines  are  0,  ±0.02,  ±0.1,  ±0.3.  The  contour  values  of  the  vorticity  axe 
±0.4,  ±0.8,  ±1.4,  ±2,  with  the  highest  magnitude  at  the  sphere  surface.  The  solid  and 
dotted  lines  in  the  figures  represent  respectively  positive  and  negative  values. 

The  pseudo-streamlines  shown  in  figures  6(a)-(f)  resemble  closely  those  for  the  inter¬ 
action  between  a  single  vortex  tube  and  a  sphere  which  was  described  in  KES.  Since  the 
description  of  the  flow  structure  with  the  aid  of  the  streamlines  is  given  in  KES,  it  will 
not  be  repeated  here,  and  only  the  vorticity  contours  will  be  described  here. 

The  vorticity  contours  in  figures  6(a)  and  6(b)  show  that  the  vortex  tubes  not  only 
are  advected  downstream  but  also  rotate  about  each  other.  The  contour  lines  of  vorticity 
in  the  figures  also  show  that  viscous  diffusion  takes  place.  It  is  well  known  that  two  co¬ 
rotating  point  vortices  located  a  distance  apart  in  an  inviscid  flow  rotate  with  constant 
angular  velocity  about  the  point  located  at  the  center  of  the  segment  connecting  them 
while  the  separation  distance  held  fixed.  On  the  other  hand,  when  two  co-rotating  vortex 
tubes  are  located  a  distance  apart  in  an  inviscid  flow  and  the  separation  distance  is  small 
enough,  they  interweave  as  well  as  rotate  about  each  other  (Zabusky,  Hughes  &  Roberts 
(1979);  Overman  and  Zabusky  (1982);  Rangel  &  Sirignano  (1989)). 
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Figures  6(c)  and  6(d)  show  that  the  vortex  tubes  contact  the  boundary  layer  of  the 
sphere  and  go  around  the  bottom  of  the  sphere.  The  reason  for  the  passage  of  the  vortex 
tubes  around  the  bottom  of  the  sphere  rather  than  around  the  top  is  as  follows.  When 
the  vortex  tubes  rotating  counter-clockwise  come  close  to  the  sphere  boundary  layer,  they 
augment  the  magnitude  of  the  vorticity  in  the  lower  boundary  layer  and  reduce  that  of  the 
vorticity  in  the  upper  boundary  layer.  Consequently,  the  vorticity  in  the  lower  boundary 
layer  induces  a  velocity  in  the  downward  direction  at  the  location  of  the  vortex  tubes  with 
higher  magnitude  than  that  induced  by  the  vorticity  in  the  upper  boundary  layer.  This 
downward  induced  velocity  advects  the  vortex  tubes  below  the  sphere  (KES). 

Figure  6(e)  shows  that  the  pairing  vortex  tubes  merge  into  one  vortex  due  to  the 
interweaving  and  the  viscosity.  Figure  6(f)  shows  that  the  vorticity  contours  around  the 
sphere  approach  that  of  the  axisymmetric  flow  as  the  tubes  are  advected  far  downstream. 

A  three-dimensional  view  of  the  pair  of  vortex  tubes  is  examined  by  considering  the  y- 
component  of  vorticity  vector.  Figures  7(a)  and  7(b)  show  the  views  of  a  three-dimensional 
contour  surface  of  =  0.2  at  t  =  6  and  t  =  21,  respectively,  for  the  flow  depicted  in 
figure  6.  The  figures  show  a  view  looking  down  with  an  acute  angle  toward  the  y-z  plane. 
The  ellipse  in  the  figures  is  the  boundary  of  the  spherical  computational  domain  viewed 
at  an  angle.  It  appears  as  a  circle  when  viewed  normal  to  the  principal  plane.  The  sphere 
is  at  the  center  of  the  domain  in  figures  7(a)  and  7(b).  Figure  7(a)  shows  that  the  two 
vortex  tubes  rotate  about  each  other.  Figure  7(b)  demonstrates  that  the  pair  of  vortex 
tubes  merge  after  some  time. 

The  resemblance  of  the  streamline  pattern  between  the  case  of  a  pair  of  vortex  tubes 
and  the  case  of  a  single  vortex  tube  indicates  that  the  force  and  moment  on  the  sphere 
due  to  a  pair  of  vortex  tubes  may  be  close  in  value  to  those  due  to  a  single  vortex  tube. 
In  the  next  subsections,  we  discuss  the  lift,  moment,  and  drag  coefficients  for  the  pair  of 
vortex  tubes  and  compare  them  with  those  for  a  single  vortex  tube. 


3.3.2  Lift,  moment,  and  drag  coefficients  and  effect  of  tube  circulation 

Figure  8  shows  the  lift  coefficients  of  the  sphere  as  a  function  of  time  for  Re  =  100, 
d0fj  =  ±1.5,  and  a  =  1.  The  lift  coefficients  are  computed  for  four  different  total 
maximum  induced  velocities  vmaxt  due  to  the  pair  of  vortex  tubes,  with  magnitudes  equal 
to  0.148,  0.295,  0.443,  and  0.590  ( vmax  =  0.1,  0.2,  0.3,  and  0.4)  normalized  by  free  stream 
velocity.  Due  to  the  sudden  placement  of  the  sphere  into  the  stream,  it  takes  a  short  time 
(0  <  t  <  0.8)  for  the  initial  flow  perturbations  to  vanish. 

When  the  pair  of  vortex  tubes  approach  the  sphere  (0  <  t  <  9),  they  produce  upwash 
resulting  in  a  positive  lift  force  on  the  sphere.  The  maximum  positive  lift  coefficient 
Cl, max i  occurs  at  about  t  =  6.8.  On  the  other  hand,  when  the  vortex  tubes  pass  the 
sphere,  they  produce  downwash  and  higher  fluid  velocity  near  the  bottom  of  the  sphere 
than  the  top  due  to  the  shear  flow  imposed  by  the  vortex  tubes  resulting  in  a  negative  lift 
force.  The  magnitude  of  the  negative  lift  is  greater  than  the  positive  lift.  The  maximum 
negative  lift  coefficient  Cl, max 2  occurs  at  about  t  =  12.2.  Cl, max  1  and  Cl, max 2  are  linearly 
proportional  to  the  total  maximum  induced  velocity.  The  maximum  positive  lift  coefficient 
Cl, maxi  is  expressed  by 

Cl, maxi  =  C  Vmaxt  i  (4) 

where  the  proportionality  constant  c  =  0.88.  The  maximum  negative  lift  coefficient 
CL,max2  is  also  expressed  by  equation  (4)  but  with  c  =  -1.62.  After  the  lift  coefficient 
reaches  its  maximum  negative  value,  it  decays  quickly  towards  zero  because  the  vortex 
tube  vorticity  is  diffused  in  the  sphere  wake.  The  time  averaged  lift  coefficient  (averaged 
over  a  time  span  between  t  =  0.8  and  the  maximum  time  24.5)  for  all  values  of  vmaxt 
is  negative  and  small  (O(10-2)).  As  mentioned  earlier,  the  behavior  of  Cx(f)  during  the 
period  0  <  t  <  0.8  is  influenced  by  the  initial  flow  perturbation,  and  thus  its  value  during 
this  initial  period  is  excluded  from  the  averaging  process.  The  root  mean  square  Cx,rm3  of 
the  lift  coefficient  as  a  function  of  time  is  also  linearly  proportional  to  vmaxt  with  c  =  0.7. 
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The  lift  coefficient  of  the  sphere  interacting  with  a  single  vortex  tube  as  a  function 
of  time  is  also  shown  as  a  reference  (marked  with  an  asterisk)  in  figure  8  for  Re  =  100, 
d0j}  -  0,  and  a  =  1  with  vmaxt  =  vmax  =  0.148.  Figure  8  shows  that  the  lift  coefficient 
of  the  sphere  interacting  with  a  pair  of  like-rotation  vortex  tubes  as  a  function  of  time 
is  approximately  the  same  as  that  of  the  sphere  interacting  with  a  single  vortex  tube  for 
the  parameters  given  above  if  the  same  total  maximum  induced  velocity  is  used  in  both 
cases.  The  dependency  of  this  phenomenon  on  the  parameters  (d0jj,  a ,  and  Re)  will  be 
discussed  in  the  following  subsections. 

Figure  9  shows  the  temporal  development  of  the  moment  coefficients  for  the  sphere 
under  the  same  conditions  as  figure  8. 

When  the  vortex  tubes  pass  the  sphere,  the  front  stagnation  point  on  the  sphere  is 
shifted  above  the  plane  x  =  0  due  to  the  downwash.  This  causes  higher  sheax  stress  in 
the  lower  left  region  compared  to  the  upper  left  region  resulting  in  a  positive  (counter¬ 
clockwise)  torque  on  the  sphere.  The  upward  shift  of  the  front  stagnation  point  also  causes 
the  shear  stress  to  be  higher  in  the  top  and  upper  right  regions  than  in  the  bottom  and 
lower  right  regions  resulting  in  a  negative  torque  on  the  sphere.  However,  the  effect  of  this 
negative  torque  is  diminished  by  the  shear  flow  induced  by  the  vortex  tubes  across  the 
sphere  which  produces  high  shear  stress  at  the  bottom  of  the  sphere.  As  a  consequence, 
a  net  high  positive  torque  acts  on  the  sphere.  The  maximum  positive  moment  coefficient 
Cm, max  occurs  at  t  =  11.5.  Cm, max  is  approximately  linearly  proportional  to  vmaxt  with 
c  =  0.11. 

When  the  vortex  tubes  approach  the  sphere  or  are  relatively  far  away  from  the  sphere, 
the  effect  of  the  shear  flow  induced  by  the  vortex  tubes  across  the  sphere  is  small,  resulting 
in  a  net  weak  torque  on  the  sphere. 

The  time  averaged  moment  coefficient  for  all  values  of  vmaxt  is  positive  and  small 
O(10-3).  The  rms  moment  coefficient  Cm,™  is  approximately  linearly  proportional  to 
Vmaxt  with  c  =  0.043.  We  note  that  the  torque  depends  only  the  distribution  of  the  shear 
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stresses  ( ttq  and  rr^)  and  is  relatively  small  compared  to  the  lift  force. 

The  moment  coefficient  of  the  sphere  interacting  with  a  single  vortex  tube  as  a  function 
of  time  is  also  shown  as  a  reference  (marked  with  an  asterisk)  in  figure  9  for  Re  =  100, 
d0jf  —  0,  and  cr  —  1  with  vmaxt  —  vmax  —  0.148.  The  pattern  of  the  moment  coefficient 
of  the  sphere  interacting  with  a  pair  of  vortex  tubes  as  a  function  of  time  is  similar  to 
that  of  the  sphere  interacting  with  a  single  vortex  tube  for  the  parameters  given  above, 
but  the  maximum  moment  coefficient  of  the  former  is  lower  than  that  of  the  latter.  This 
shows  that  the  moment  coefficient  is  more  sensitive  to  the  offset  distance  than  the  lift 
coefficient.  This  will  be  discussed  in  the  next  section  in  detail. 

Figure  10  shows  the  drag  coefficients  of  the  sphere  as  a  function  of  time  for  the  same 
conditions  as  figure  8.  The  drag  coefficients  are  computed  for  four  different  values  of 
v„„Tt  as  in  figure  8,  in  addition  to  vmaxt  =  0  which  corresponds  to  the  axisymmetric  flow 
without  the  vortex  tubes. 

As  discussed  earlier,  the  sudden  placement  of  the  sphere  in  the  flow  results  in  initially 
large  values  of  shear  stress  and  pressure  on  the  sphere,  and  hence  a  large  drag  as  shown  in 
figure  10.  When  the  vortex  tubes  approach  the  sphere,  the  pressure  at  the  front  stagnation 
point  is  lower  than  that  of  the  axisymmetric  flow  past  a  sphere  due  to  the  low  pressure 
at  the  center  of  the  vortex  tube.  Also,  the  maximum  shear  stresses  in  the  upper  and 
lower  regions  of  the  sphere  are  lower  than  those  of  the  axisymmetric  flow.  This  causes 
the  drag  on  the  sphere  to  be  lower  than  that  of  the  axisymmetric  flow  without  the  vortex 
tube.  As  the  vortex  tubes  move  around  the  bottom  of  the  sphere,  the  front  stagnation 
point  is  shifted  above  the  plane  x  =  0  due  to  the  downwash.  Consequently,  high  pressure 
and  high  shear  stress  act  in  the  upper  and  lower  left  regions,  respectively.  This  increases 
the  drag  during  the  period  9  <  t  <  13.4.  For  t  >  13.4,  the  drag  approaches  that  of 
the  axisymmetric  flow  as  the  vortex  tube  moves  further  downstream.  The  time  averaged 
value  of  the  deviation  of  the  drag  coefficient  from  that  of  the  axisymmetric  flow  past  a 
sphere  for  all  values  of  vmax  is  nearly  zero  (O(10-4)). 
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The  pattern  of  the  drag  coefficient  of  the  sphere  interacting  with  a  pair  of  vortex 
tubes  as  a  function  of  time  is  similar  to  that  of  the  sphere  interacting  with  a  single  vortex 
tube  for  the  parameters  given  above  (compared  with  figure  11  in  KES),  but  the  largest 
deviation  of  the  former  from  the  case  of  the  axisymmetric  flow  occurs  earlier  than  that 
of  the  latter.  The  reason  is  that  due  to  the  rotation  about  each  other,  one  of  the  vortex 
tubes  in  the  former  approaches  the  sphere  faster  than  the  vortex  tube  in  the  latter. 

3.3.3  Effects  of  the  size  and  the  offset  distance  of  the  vortex  tubes 

The  effects  of  the  size  of  the  vortex  tubes  on  the  flow  field  are  studied  by  performing 
computations  similar  to  those  in  the  previous  section  for  Re  =  100,  d0fj  =  ±1.5,  and  five 
different  sizes  of  the  vortex  tubes,  a  =  0.25,  0.5,  2,  3,  and  4  in  addition  to  the  base  case 
<7=1. 

Table  1  shows  CL,m „u  CL,max2 ,  CL,rm„  Cm, max,  and  CMrms  as  a  function  of  the  vortex 
tube  size  which  covers  six  different  initial  radii  of  the  vortex  tube,  a  =  4,3,2, 1,0.5,  and 
0.25,  for  vma xt  =  0.1.  Another  computation  with  different  vmaxi  showed  that  all  the  lift 
and  moment  coefficients  are  linearly  proportional  to  vmaxt  at  each  a.  When  a  >  2,  Cl, maxi 
and  CL,rm»  become  independent  of  <r,  but  the  magnitudes  of  Ci„max2,  Cm, max,  and  C\f,rms 
for  <j  =  4  axe  smaller  than  those  for  <r  =  2  and  3.  When  a  approaches  zero,  all  the 
coefficients  tend  to  be  proportional  to  (a  vma£t)  or  (cvmax)  which  is  proportional  to  the 
circulation  of  the  vortex  tube.  For  example,  Cl,™*  is  expressed  by 

Cl  ,rm«  —  ^1  ^ maxt  ?  2  ^  (T  ^  4 

=  C2  vmaxt  <7n,  0.25  <  <7  <  2,  0.75  >  n  >  0.3  ,  (5) 

where  the  constant  C\  =  1  and  c2  =  0.7,  and  n  depends  on  o  and  should  approach  unity 
as  reaches  zero.  For  Cl, maxi ,  1.1  and  c2  0.38.  ^l,thch2,  C*7^,max>  and  CM,Tm8  for 

<7  <  3  axe  also  expressed  by  equation  (5)  with  c\  =  —2  and  c2  =  —1.65,  cj  =  0.13  and 
c2  =  0.11,  and  C\  —  0.053  and  c2  =  0.04,  respectively.  The  time  averaged  value  of  the 
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deviation  of  the  drag  coefficient  from  that  of  the  axisymmetric  flow  past  a  sphere  for  all 
values  of  a  is  nearly  zero  (O(10-4)). 

Comparing  the  results  in  table  1  (based  on  vmaxt )  with  those  of  the  sphere  interacting 
with  a  single  vortex  tube  (KES,  table  4  (based  on  umax)),  it  is  found  that  the  magnitudes 
of  the  lift  coefficients  in  table  1  are  within  2%  to  20%  of  those  in  KES,  with  the  largest 
deviation  occurring  at  a  =  0.25. 

Note  that  Cl, maxi,  Cm, maxi  and  Cm ,rm»  for  a  =  4  axe,  respectively,  smaller  than  those 
for  cr  =  2  and  3  due  to  the  shear  flow  effect  explained  in  KES. 

Now,  the  effects  of  the  offset  distance  on  the  flow  field  are  investigated  by  varying  d0jj 
for  Re  =  100  and  a  =  4.  The  computation  was  performed  for  d0jj  =  0,  ±1,  ±2,  ±3,  and 
±4  in  addition  to  the  base  case  d0jj  =  ±1.5.  Note  that  the  case  of  d0/f  =  0  corresponds 
to  the  interaction  between  a  single  vortex  tube  and  a  sphere. 

It  is  found  that  Cl, maxi ,  C L,maxh  CL,rmn  Cm, maxi  and  C Mrma  for  each  dQf  f  are  linearly 
proportional  to  vmaxt  as  in  the  case  of  d0jj  =  ±1.5.  The  triangular  symbols  in  figure  11 
show  CL,rms  as  a  function  of  \d0jj\  for  Re  =  100  and  a  =  4  while  the  maximum  induced 
velocity  (or  the  circulation)  of  each  vortex  tube  is  kept  as  a  constant,  vmax  =  0.2.  The 
triangular  symbols  show  that  C^rms  decays  rapidly  as  >  0.  On  the  other  hand, 

the  circular  symbols  in  figure  11  show  Ci„rm,  as  a  function  of  |d0//|  for  Re  =  100  and 
<t  =  4  while  the  total  maximum  induced  velocity  due  to  the  two  vortex  tubes  is  kept  as 
a  constant,  vmaxt  =  0.2.  The  circular  symbols  show  that  the  magnitudes  of  the  rms  lift 
coefficients  for  dQ //  =  ±1,  ±1.5,  ±2,  ±3,  and  ±4  axe  close  to  that  for  d0jj  =  0.  The 
behavior  of  Cl, maxi  and  Cl, max 2  as  a  function  of  |<f0//|  is  similar  to  that  of  Cl ,rm«- 

Examination  of  the  effect  of  the  offset  distance  for  <7  =  1  and  2  shows  that  the  lift 
coefficient  of  the  sphere  interacting  with  a  pair  of  vortex  tubes  as  a  function  of  time  is 
nearly  identical  to  that  of  the  sphere  interacting  with  a  single  vortex  tube  if  the  separation 
distance  between  the  tube  centers  is  less  than  2  i/a  vortex  tube  diameter  for  Re  =  100 
and  vmaxt  instead  of  vmax  is  used  in  the  former. 
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The  triangular  symbols  and  the  circular  symbols  in  figure  12  show  Cm,™*  as  a  function 
of  \d0fj\  for  the  same  parameters  as  used  for  Cx,rm4.  Figure  12  shows  that  the  magnitude 
of  the  rms  moment  coefficient  decays  more  rapidly  than  that  of  the  rms  lift  coefficient  as 
d0j  <  increases.  The  behavior  of  Cm, max  as  a  function  of  dafj  is  similar  to  that  of  Cm, rms- 
The  moment  coefficient  of  the  sphere  interacting  with  a  pair  of  vortex  tubes  as  a  function 
of  time  is  nearly  identical  to  that  of  the  sphere  interacting  with  a  single  vortex  tube  if 
the  separation  distance  between  the  tube  centers  is  less  than  y/a  vortex  tube  diameter 
for  Re  =  100  and  vmaxt  instead  of  vmax  is  used  in  the  former. 

3.3.4  Effects  of  Reynolds  number 

Computations  like  those  in  section  3.3.2  are  made  for  four  different  Reynolds  numbers 
in  the  range  of  20  <  Re  <  80,  d0fj  =  ±1.5,  and  1  <  a  <  4  in  addition  to  the  base  case 
Re  =  100. 

A  result  like  that  shown  in  section  3.3.3  for  Re  =  100  is  obtained.  Cl, maxi  and  CL,rm$ 
axe  linearly  proportional  only  to  vmaxt  and  independent  of  a  when  c  >  2  at  fixed  Reynolds 
number  as  in  the  case  of  Re  =  100.  Cl, maxi  dependence  on  Reynolds  number  may  be 
expressed  by 

Cl, maxi  —  A  Vmaxt  die  ,  (6) 

where  A  =  8.9  and  P  =  —0.45  for  2  <  a  <  4.  Cl,™  may  be  also  expressed  by  equation  (6) 
with  A  =  8.1  and  P  =  —0.45  for  2  <  cr  <  4.  Cm, max  and  Cm,™  may  be  also  represented 
by  equation  (6)  with  A  =  5.5  and  P  —  —0.83  for  the  former,  and  A  —  3.1  and  P  =  —0.88 
for  the  latter  for  2  <  a  <  3. 

Now,  the  effect  of  the  offset  distance  for  20  <  Re  <  80  in  addition  to  the  base  case 
Re  =  100  is  discussed. 

The  triangular  symbols  in  figure  13  show  Cl, rms  as  a  function  of  \d0jj\  for  Re  =  20 
and  cr  =  4  while  the  maximum  induced  velocity  (or  the  circulation)  of  each  vortex  tube  is 


20 


kept  as  a  constant,  vmax  —  0.2.  The  triangular  symbols  show  that  C^rms  decays  rapidly 
as  \d0ff  \  >  0.  On  the  other  hand,  the  circular  symbols  show  Cl, rms  as  a  function  of  \d0ff\ 
for  Re  =  20  and  a  =  4  while  the  total  maximum  induced  velocity  due  to  the  two  vortex 
tubes  is  kept  as  a  constant,  vmaxt  =  0.2.  The  circular  symbols  show  that  the  magnitudes 
of  the  rms  lift  coefficients  for  dQjj  =  ±2  and  ±4  are  close  to  that  for  d0/j  =  0.  The 
behavior  of  Cl, maxi  and  Ci,max2  as  a  function  of  \d0jj\  resembles  that  of  CL,rms- 

The  results  for  the  range  of  a  values  indicate  that  the  lift  coefficient  of  the  sphere 
interacting  with  a  pair  of  like-rotation  vortex  tubes  as  a  function  of  time  is  nearly  identical 
to  that  of  the  sphere  interacting  with  a  single  vortex  tube  if  the  separation  distance 
between  the  tube  centers  is  less  than  2  yfa  vortex  tube  diameter  for  Re  =  20  and  vmaxt 
instead  of  vmax  is  used  in  the  former  case.  The  same  result  as  above  was  obtained  at 
different  Reynolds  numbers,  Re  =  40,  60,  and  80. 

The  triangular  symbols  and  the  circular  symbols  in  figure  14  show  CM.rms  as  a  function 
of  \d0fj\  for  the  same  parameters  as  used  for  Cl, rms-  The  figure  shows  the  magnitude  of 
the  rms  moment  coefficient  decays  more  rapidly  than  that  of  the  rms  lift  coefficient  as 
d0jf  increases.  The  behavior  of  Cm, max  as  a  function  of  daff  resembles  that  of  Cm, rms-  It 
is  found  that  the  moment  coefficient  of  the  sphere  interacting  with  a  pair  of  vortex  tubes 
as  a  function  of  time  is  nearly  identical  to  that  of  the  sphere  interacting  with  a  single 
vortex  tube  if  the  separation  distance  between  the  tube  centers  is  less  than  yfa  vortex 
tube  diameter  for  Re  —  20  and  vmaxt  instead  of  vmax  is  used  in  the  former  case.  The  same 
result  as  above  was  obtained  at  different  Reynolds  numbers,  Re  =  40,  60,  and  80. 

In  summary,  the  comparison  of  the  results  from  this  section  with  those  from  the 
previous  section  shows  that  the  range  of  the  offset  distance  for  which  the  lift  and  moment 
coefficients  of  the  sphere  interacting  with  a  pair  of  vortex  tubes  are  nearly  identical  to 
those  of  the  sphere  interacting  with  a  single  vortex  tube  is  independent  of  Reynolds 
number  for  20  <  Re  <  100. 
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3.4  Interactions  of  an  array  of  vortex  tubes  of  like  rotation  and 
a  moving  sphere 

Our  results  for  the  cases  of  a  fixed  spherical  particle  interacting  with  a  single  advecting 
vortex  and  with  an  advecting  pair  of  vortices  can  now  be  used  to  calculate  the  trajectory 
of  a  moving  spherical  particle  interacting  with  an  array  of  like-rotating  vortex  tubes. 
Opposite-rotating  vortices  are  less  interesting  since  they  produce  no  lift  or  deflection. 
Studying  these  interactions  can  improve  our  understanding  of  the  behavior  of  a  particle 
(or  droplet)  interacting  with  eddies  of  comparable  length  scale  in  a  turbulent  flow.  For 
this  study,  the  array  will  be  a  linear  arrangement  of  single  vortices.  Since  a  single  vortex 
and  a  pair  of  like-rotating  vortices  produce  comparable  effects  on  the  sphere,  this  choice 
should  not  be  critical. 

Figure  15  shows  the  initial  flow  geometry  where  a  spherical  particle  is  injected  into  an 
array  of  infinite  number  of  counter-clockwise  rotating  vortices  which  are  located  on  the 
negative  z-axis  with  center-to-certer  nondimensional  distance  of  24.  Since  the  life  time  of 
a  vortex  tube  is  short  (2i rcr /vmax)  compared  to  the  travel  time  (or  life  time)  of  the  particle 
(or  droplet),  it  is  assumed  that  the  next  vortex  with  the  same  strength  as  the  first  vortex 
is  generated  when  the  sphere  passes  the  first  vortex.  The  deflection  of  the  moving  sphere 
will  cause  the  offset  distance  to  vary  from  one  collision  with  a  vortex  tube  to  the  next. 
Therefore,  d0fj  is  a  time-dependent  quantity. 

We  assume  that  the  particle  is  constrained  to  move  only  in  the  x-z  plane.  The  aim 
is  to  calculate  the  trajectory  in  the  two-dimensional  plane  from  the  already-known  time 
evolution  of  Cc(f)  and  Cx(t)  as  given  by  KES.  The  particle  trajectory  as  a  function  of 
time  is  computed  by  solving  the  following  system  of  two  ordinary  differential  equations 
which  are  the  nondimensional  form  of  the  Newton’s  equation  of  motion  in  the  z  and  x 
directions. 
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(7a) 


=  T~(~cD(t)  c»3»  -  cut)  sine)  (Ul  +  ul) 

at  o  pr 

=  g ^-(-C'u(i)  sinO  +  CL(t)  cos9)  (U2  +  U2) ,  (7b) 

where  Uz  and  Ux  are,  respectively,  the  sphere  velocities  in  z  and  x  direction,  tan6  = 
Ux/{—Uz),  and  pT  is  the  ratio  of  the  particle  density  to  the  fluid  density.  Initially,  0  =  0. 
The  term  {U^+U2)  arises  since  the  velocities  are  normalized  by  the  initial  particle  velocity 
while  the  drag  and  lift  coefficients  are  normalized  by  the  instantaneous  particle  velocity. 
The  gravity  force  is  neglected  in  this  formulation. 

Our  numerical  results  (KES)  for  the  time-dependent  lift  coefficient  of  a  spherical  par¬ 
ticle  interacting  with  a  single  vortex  tube  depends  on  the  offset  distance  d0fj ,  the  vortex 
core  size  cr,  and  Reynolds  number  Re  and  can  be  summarized  as: 

—  A  (Ct(f)]t,  exp(^~)  f(Re,vmax)  for  C  ^  <  d0jt  <  D  (8) 

\C 

where  the  following  combinations  of  values  apply: 


A 

B 

C 

D 

1.15 

0.3 

—  CO 

-0.7 

1 

0.1 

-0.7 

0 

1 

0 

0 

1 

e0-1 

-0.1 

1 

1.7 

1.15  e0-3 

0.3 

1.7 

oo 

Also,  [ Ci{t)]b  is  the  lift  coefficient  for  the  base  case  where  Re  =  100  and  d0jj  =  0  with 
Vmax  =  vmax,b,  f(R^A ’max)  =  (100/i?e)m  vmax/ vmaXtb,  and  the  exponent  m  is  given  by 
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m 

a 

t 

0.375 

i 

0  <  t  <  9 

0.45 

i 

9  <  t  <  24 

0.44 

2  <<r<4 

0  <t  <  9 

0.51 

2  <  <x  <  4 

9  <  t  <  24 

The  time-averaged  drag  coefficient  of  the  sphere  in  the  flow  with  a  vortex  tube  differs 
by  0.01%  to  5%  (depending  on  the  offset  distance)  from  the  time- averaged  drag  coefficient 
of  the  axisymmetric  flow  for  Re  =  100  and  vmax  =  0.2.  (Refer  to  equation  (19)  in  KES.) 
Therefore,  the  time-dependent  drag  coefficient  in  equations  (7a)  and  (7b)  is  approximated 
by  the  time-dependent  drag  coefficient  obtained  from  the  axisymmetric  flow  generated  by 
a  spherical  particle  injected  into  a  quiescent  fluid.  The  spherical  particle  in  this  flow 
experiences  the  drag  force  and  thus  is  retarded.  This  axisymmetric  drag  coefficient  was 
computed  as  a  function  of  time  and  instantaneous  Reynolds  number  by  using  the  code 
which  has  been  developed  for  the  time-dependent  axisymmetric  flow.  The  torque  on  the 
sphere  is  neglected  since  the  magnitude  of  the  moment  coefficient  is  small  and  less  than 
8%  of  the  lift  coefficient  magnitude  (table  4  of  KES). 

Figure  16  shows  two  trajectories  of  the  sphere  during  the  dimensionless  time  period 
between  0  and  24  for  initial  particle  Reynolds  numbers  50  and  100  with  density  ratio  200 
(which  is  the  ratio,  for  example,  of  n-octane  density  to  that  of  air  under  10  atmospheres 
of  pressure).  The  initial  vortex  size  is  three  times  the  sphere  radius,  and  the  initial  offset 
distance  of  the  sphere  is  zero.  The  initial  maximum  induced  velocity  of  the  vortex  tube  is 
0.2  normalized  by  the  initial  sphere  velocity.  The  sphere  initially  moves  upward  due  to  the 
vortex  up  wash  and  then  moves  downward  due  to  the  vortex  downwash.  The  maximum 
positive  deflection  for  the  case  of  Re0  =  50  is  higher  than  that  of  Re0  =  100. 

Figure  17  shows  two  trajectories  of  the  sphere  which  are  traced  from  the  initial  injection 
beyond  the  time  period  of  figure  16  until  particle  Reynolds  number  reaches  unity  for 
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initial  particle  Reynolds  numbers  50  and  100  with  the  same  initial  parameters  as  in 
figure  16.  Since  the  final  Reynolds  number  is  small,  these  trajectories  are  approximately 
those  corresponding  to  the  whole  particle  motion  until  it  stops  relative  to  the  fluid.  A 
counter-clockwise  rotating  vortex  tube  produces  not  only  upwash  downstream  of  itself 
and  downwash  upstream  of  itself  but  it  also  causes  a  shear  flow  across  the  sphere  when 
it  passes  the  sphere.  The  combined  effect  of  the  downwash  and  the  shear  flow  causes  the 
magnitude  of  the  maximum  negative  lift  to  be  greater  than  the  maximum  positive  lift 
magnitude.  Therefore,  the  average  lift  coefficient  averaged  over  the  time  span  24  (the 
interaction  time  with  one  vortex  tube)  is  one  order  of  magnitude  less  than  the  rms  lift 
coefficient  and  negative  due  to  the  shear  flow  effect.  This  small  negative  value  of  the 
average  lift  coefficient  becomes  important  when  the  sphere  interacts  with  an  array  of 
many  vortices.  Thus,  the  sphere  travels  upward  only  for  the  short  initial  time  period  and 
then  moves  downward  for  the  most  of  the  time  until  it  stops.  The  final  deflection  ratios 
defined  by  the  ratio  of  the  final  position  xj  to  zj  of  the  sphere  are  1/36  for  the  case  of 
Re0  =  100  and  1/34  for  the  case  of  Re0  =  50.  However,  the  final  deflection  for  the  case  of 
Re0  =  100  is  higher  than  that  of  Re0  =  50,  because  the  sphere  for  the  case  of  Re0  =  100 
possesses  higher  initial  momentum  and  it  travels  farther  than  that  of  Re0  =  50. 

Figure  18  shows  four  trajectories  of  the  sphere  during  the  dimensionless  time  period 
between  0  and  24  for  the  density  ratio  25,  50,  100,  and  200  with  Reynolds  number  100 
and  the  same  parameters  for  <7,  d0ff,  and  vmax  as  in  figure  16.  A  sphere  with  lower  density 
ratio  initially  deflects  more  than  a  sphere  with  higher  density  ratio  as  shown  in  figure  19. 
Four  trajectories  of  the  sphere  which  are  traced  from  the  initial  injection  until  particle 
Reynolds  number  reaches  unity  for  the  density  ratio  25,  50,  100,  and  200  with  the  same 
initial  parameters  as  in  figure  18.  However,  the  final  transverse  displacement  increases 
with  density  ratio  because  the  sphere  with  higher  density  ratio  possesses  higher  initial 
momentum  and  it  travels  farther  than  the  sphere  with  lower  density  ratio. 

The  larger  vmax  causes  the  larger  sphere  deflection;  however,  the  sphere  deflection  is 
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not  linearly  proportional  to  vmax  due  to  the  nonlinearity  of  the  equations  (7a)  and  (7b). 

The  results  of  figures  16  and  18  indicate  that  the  sphere  would  experience  slightly 
lower  drag  than  that  of  a  sphere  subjected  to  an  axisymmetric  flow  when  it  passes  the 
first  vortex  tube.  This  lower  drag  is  caused  by  the  upward  motion  of  the  sphere  due 
to  the  upwash  of  the  approaching  vortex  tube,  and  thus  the  center  of  the  vortex  tube  is 
located  below  the  front  stagnation  point  of  the  sphere.  This  causes  lower  dynamic  pressure 
ahead  of  the  front  stagnation  point.  However,  the  sphere  would  experience  higher  drag 
eventually  when  it  passes  more  than  one  vortex  tube  and  travels  downward.  Due  to  the 
downward  motion,  the  vortex  tubes  are  located  above  the  front  stagnation  point  of  the 
sphere,  causing  higher  dynamic  pressure  ahead  of  the  front  stagnation  point. 

The  original  computations  in  KES  and  in  section  3.3  were  made  for  a  nondimensional 
time  duration  of  24  and  24.5,  respectively.  However,  the  trajectory  calculations  presented 
in  figures  17  and  19  use  the  basic  information  from  those  original  computations  to  yield 
trajectory  predictions  for  much  longer  periods. 

4  Conclusions 

In  order  tc  improve  +he  understanding  of  the  physics  of  interaction  between  a  particle 
and  eddies  of  comparab1  -  length  scale  in  a  carrier  flow,  the  unsteady,  three-dimensional, 
incompressible,  viscous  flow  interactions  between  a  pair  of  vortex  tubes  advected  by  a 
uniform  free  stream  and  a  spherical  particle  suddenly  placed  and  held  fixed  in  space  were 
investigated  numerically  for  a  range  of  particle  Reynolds  number  20  <  Re  <  100. 

When  the  top  and  bottom  vortex  tubes  have  positive  and  negative  circulations,  re¬ 
spectively,  the  magnitude  of  the  induced  velocity  due  to  the  vortex  tubes  is  added  to 
the  base  flow  velocity  along  the  stagnation  streamline.  This  causes  the  pressure  at  the 
stagnation  point  and  the  shear  stresses  in  the  upper  and  lower  left  regions  to  be  higher 
than  those  of  the  axisymmetric  flow  past  a  sphere,  thus  increasing  the  drag.  On  the  other 
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hand,  when  the  top  and  bottom  vortex  tubes  have  negative  and  positive  circulations, 
respectively,  the  induced  velocity  due  to  the  vortex  tubes  is  subtracted  from  the  base  flow 
velocity  along  the  stagnation  streamline.  This  causes  the  pressure  at  the  stagnation  point 
and  the  shear  stresses  in  the  upper  and  lower  left  regions  to  be  lower  than  those  of  the 
axisymmetric  flow  past  a  sphere,  thus  reducing  the  drag.  The  lift  and  moment  axe  zero 
for  this  symmetric  configuration. 

The  interactions  between  a  sphere  and  like-rotating  a  pair  of  cylindrical  vortex  tubes 
initially  located  ten  radii  upstream  from  the  center  of  the  sphere  were  investigated.  The 
lift  and  moment  coefficients  of  the  sphere  interacting  with  a  pair  of  vortex  tubes  as  a 
function  of  time  are  nearly  identical,  respectively,  to  those  of  the  sphere  interacting  with 
a  single  vortex  tube  if  the  separation  distance  between  the  tube  centers  is  less  than  2  yfd 
vortex  tube  diameter  for  the  lift  coefficient  and  less  than  \fc  vortex  tube  diameter  for 
the  moment  coefficient;  here,  vmaxt  instead  of  vmax  is  used  in  the  case  of  a  pair  of  vortex 
tubes,  where  vmax  is  the  maximum  induced  velocity  due  to  one  vortex  without  presence  of 
the  other  and  vmaxt  is  the  total  maximum  induced  velocity  due  to  the  pair  of  vortices.  In 
particular,  lift  and  moment  coefficients  axe  linearly  proportional  to  the  maximum  induced 
velocity.  The  moment  coefficient  is  negligible  compared  to  the  lift  coefficient. 

The  two-dimensional  trajectories  of  a  spherical  particle  interacting  with  an  array  of 
vortices  whose  sizes  are  comparable  to  the  sphere  size  have  been  examined.  The  time- 
dependent  drag  and  lift  forces  (KES)  for  the  case  of  a  spherical  particle  interacting  with  a 
single  vortex  were  used  to  calculate  the  two-dimensional  trajectory  of  a  moving  spherical 
particle  interacting  with  an  array  of  vortex  tubes  of  like  rotation.  The  present  results 
show  that  the  shear  flow  across  the  sphere  induced  by  a  vortex  tube  is  responsible  for 
the  net  deflection  of  a  sphere  interacting  with  an  array  of  vortex  tubes.  Thus,  the  sphere 
eventually  deflects  in  the  direction  of  increasing  relative  velocity.  The  deflection  ratio 
(ratio  of  sphere  final  location  in  z  and  x  directions)  of  the  sphere  increases  with  decreasing 
initial  Reynolds  number  and  with  decreasing  density  ratio.  However,  the  total  deflection 
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increases  with  increasing  the  initial  Reynolds  number  and  the  density  ratio  because  higher 
momentum  causes  the  sphere  to  travel  farther. 

A  turbulent  flow  possesses  a  wide  spectrum  of  eddy  sizes.  In  order  to  enhance  the 
understanding  of  particle  motion  in  a  turbulent  flow,  we  are  currently  investigating  the 
motion  of  a  freely  moving  particle  interacting  with  a  large  vortex  tube  whose  size  is  of 
the  order  of  an  integral  length  scale,  i.e.  at  the  other  end  of  the  spectrum  relative  to  the 
present  case.  Preliminary  results  of  the  large  vortex  study  show  that  the  lift  force  on  the 
particle  is  much  smaller  than  that  due  to  the  small  vortices  (in  the  present  paper)  with 
the  same  maximum  induced  velocity  (Kim,  Elghobashi  &  Sirignano  (1996)).  Our  direct 
solution  of  the  three-dimensional  Navier-Stokes  equations  over  a  freely  moving  particle 
confirms  that  due  mainly  to  the  drag  force,  the  particle  travels  in  a  curved  trajectory 
that  depends  on  the  direction  of  vortex  rotation  and  the  Stokes  number  of  the  particle. 
This  indicates  that  the  mechanism  of  particle  dispersion  due  to  the  interaction  with  small 
vortices  is  quite  different  from  that  due  to  the  interaction  with  a  large  vortex. 

Acknowledgement 

This  work  has  been  supported  by  the  Air  Force  Office  of  Scientific  Research  under 
grant  No.  F49620-93-1-0028  with  Dr.  Julian  Tishkoff  acting  as  the  technical  monitor.  We 
would  like  to  thank  Mr.  Lyle  Wiedeman  for  his  assistance  in  using  a  three-dimensional 
graphic  package  Application  Visualization  System  (AVS).  The  support  of  the  San  Diego 
Supercomputer  Center  and  the  San  Diego  Supercomputer  Center  under  a  block  grant  of 
the  Office  of  Academic  Computing  of  UCI  are  gratefully  appreciated. 


28 


5  References 


Anderson,  D.  A.,  Tannehill,  J.  C.  k  Pletcher,  R.  H.  1984  Computational  Fluid  Mechanics 
and  Heat  Transfer.  Hemisphere  Publishing. 

Clift,  R.,  Grace,  J.  R.  k  Weber,  M.  E.  1978  Bubbles,  Drops,  and  Particles.  Academic 
Press,  New  York. 

Kim,  I.,  Elghobashi  S.  &  Sirignano,  W.  A.  1995  Unsteady  flow  interactions  between  an 
advected  cylindrical  vortex  tube  and  a  spherical  particle.  J.  Fluid  Mech.  288,  123-155. 

Kim,  I.,  Elghobashi  S.  k  Sirignano,  W.  A.  1996  The  motion  of  a  spherical  particle  in 
unsteady  flows  at  moderate  Reynolds  Numbers  3fth  AIAA  Aerospace  Sciences  Meeting 
Preprint  96-0081. 

Overman,  E.  A.  k  Zabusky,  N.  J.  1982  Evolution  and  merger  of  isolated  vortex  structures. 
Phys.  Fluids  25,  1297-1305. 

Patnaik,  G.  1986  A  numerical  solution  of  droplet  vaporization  with  convection.  Ph.D. 
Dissertation,  Carnegie-Mellon  University. 

Patnaik,  G.,  Sirignano,  W.  A.,  Dwyer,  H.  A.  k  Sanders,  B.  R.  1986  A  numerical  technique 
for  the  solution  of  a  vaporizing  fuel  droplet.  Prog.  Astro.  Aero.  105,  253-266. 

Rangel,  R.  k  Sirignano,  W.  A.  1989  The  dynamics  of  vortex  pairing  and  merging.  27th 
AIAA  Aerospace  Sciences  Meeting  Preprint  89-0128. 

Spalart,  P.  R.  1982  Numerical  simulation  of  separated  flows.  Ph.  D.  Thesis,  Stanford 
University. 

Taneda,  S.  1956  Experimental  investigation  of  the  wake  behind  a  sphere  at  low  Reynolds 
number.  J.  Phys.  Soc.  Japan  11  1104-1108. 


29 


Vinokur,  M  1983  On  one-dimensional  stretching  functions  for  finite- difference  calculations. 
J.  Comput.  Phys.  50,  215-234. 

Zabusky,  N.  J.,  Hughes  M.  &  Roberts  K.  V.  1979  Contour  dynamics  for  the  Euler  equa¬ 
tions  in  two  dimensions.  J.  Comp.  Phys.  30,  96-106. 


a 
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0.111 
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0.103 

0.011 

0.0051 
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0.111 
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0.102 

0.013 

0.0053 
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0.094 
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0.0053 
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0.011 

0.0043 

0.5 

0.058 

-0.116 

0.046 

0.0065 

G.0023 

0.25 

0.034 

-0.070 

0.027 

0.0034 

0.0012 

Table  1.  Maximum  positive  and  negative  lift  coefficients, 
rms  lift  coefficient,  maximum  moment  coefficient, 
and  rms  moment  coefficient  as  a  function  of 
the  size  of  vortex  tube  for  Re  =  100  and 
d0ff  =  1.5  with  vmaxi  =  0.1. 


Figure  Captions 


Figure  1.  Flow  geometry  and  coordinates 

Figure  2.  Pseudo-streamlines  (left  column)  and  contour  lines  of  y-component  vorticity 
(right  column)  in  the  principal  plane  at  (a)  t  =  0,  (b)  3,  (c)  6,  (d)  9,  (e)  12, 
and  (f)  15  for  Re  =  100,  d0fj  =  ±1.5,  <7  =  1,  and  vmaxt  =  0.738 
with  top-positive  and  bottom-negative  circulations. 

Figure  3.  Drag  coefficients  of  the  sphere  as  a  function  of  time  and  vmaxt 
for  Re  =  100,  d0jj  =  ±1.5,  and  <7=1. 
with  top-positive  and  bottom-negative  circulations. 

Figure  4.  Pseudo-streamlines  (left  column)  and  contour  lines  of  y-component  vorticity 
(right  column)  in  the  principal  plane  at  (a)  t  =  0,  (b)  3,  (c)  6,  (d)  9,  (e)  12, 
and  (f)  15  for  Re  =  100,  d0fj  =  ±1.5,  <7  =  1,  and  vmaxt  =  0.738 
with  top-negative  and  bottom-positive  circulations. 

Figure  5.  Drag  coefficients  of  the  sphere  as  a  function  of  time  and  vmaxt 
for  Re  =  100,  d0ff  =  ±1.5,  and  a  =  1. 
with  top-negative  and  bottom-positive  circulations. 

Figure  6.  Pseudo-streamlines  (left  column)  and  contour  lines  of  y-component  vorticity 
(right  column)  in  the  principal  plane  at  (a)  t  =  1,  (b)  6,  (c)  10,  (d)  15, 

(e)  21,  and  (f)  30  for  Re  —  100,  d0jj  =  ±1.5,  <7=1, 
and  n ti  —  0.59. 

Figure  7.  A  view  of  three-dimensional  contour  surfaces  of  uy  =  0.2  at  (a)  t  =  6 
and  (b)  t  =  21  for  the  flow  depicted  in  figure  2 

Figure  8.  Lift  coefficients  of  the  sphere  as  a  function  of  time  and  vmaxt 
for  Re  =  100,  <f0//  =  ±1.5,  and  <7=1. 


Figure  9.  Moment  coefficients  of  the  sphere  under  the  same  conditions  as  figure  7. 

Figure  10.  Drag  coefficients  of  the  sphere  under  the  same  conditions  as  figure  7. 

Figure  11.  Rms  lift  coefficients  of  the  sphere  as  a  function  of  \d0jj \ 
for  Re  =  100  and  a  —  4. 

Figure  12.  Rms  moment  coefficients  of  the  sphere  as  a  function  of  \d0jj\ 
for  Re  =  100  and  a  =  4. 

Figure  13.  Rms  lift  coefficients  of  the  sphere  as  a  function  of  \d0jj\ 
for  Re  =  20  and  a  —  4. 

Figure  14.  Rms  moment  coefficients  of  the  sphere  as  a  function  of  \d0/j\ 
for  Re  =  20  and  a  =  4. 

Figure  15.  Initial  flow  geometry  for  a  sphere  injected  into  an  array  of 
infinite  number  of  vortex  tubes. 

Figure  16  Two  trajectories  of  the  sphere  during  the  time  period  between  0  and  24 
for  initial  Reynolds  numbers  50  and  100  with  density  ratio  200. 

Figure  17  Two  trajectories  of  the  sphere  traced  from  the  initial  injection  until 

Reynolds  number  reaches  unity  for  initial  Reynolds  numbers  50  and  100  with 
with  density  ratio  200. 

Figure  18  Four  trajectories  of  the  sphere  during  the  time  period  between  0  and  24 

for  the  density  ratio  25,  50,  100,  and  200  with  initial  Reynolds  number  100. 

Figure  19  Four  trajectories  of  the  sphere  traced  from  the  initial  injection  until 

Reynolds  number  reaches  unity  for  the  density  ratios  25,  50,  100,  and  200 
with  initial  Reynolds  number  100. 
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The  Motion  of  A  Spherical  Particle  in  Unsteady  Flows  at 

Moderate  Reynolds  Numbers  * 
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Abstract 

The  equations  governing  the  motion  of  a  spherical 
particle  proposed  by  previous  workers  are  examined 
and  compared  with  the  results  of  the  numerical  solu¬ 
tion  of  the  full  Navier-Stokes  equations  for  unsteady, 
axisymmetric  flow  around  a  freely  moving  sphere  ini¬ 
tially  injected  into  an  oscillating  flow  and  for  unsteady, 
three-dimensional  flow  around  a  freely  moving  sphere 
interacting  with  a  large  vortex  tube.  As  a  result,  we 
propose  a  modified  equation  of  the  particle  motion  and 
demonstrate  its  superiority  to  the  previously  proposed 
equations  for  both  rectilinear  and  two-dimensional  mo¬ 
tion  over  a  wide  range  of  Reynolds  number  and  of  den¬ 
sity  ratio. 

Nomenclature 

a  sphere  radius 

Cd  drag  coefficient 

Cl  lift  coefficient 

Cm  moment  coefficient 

g  gravity  vector 

mj  mass  of  the  fluid  displaced  by  the  sphere 

rrip  mass  of  the  spherical  particle 

jVi,  X2,  N3  numbers  of  grids  in  £,  r)X  directions 

Ret  Reynolds  number  based  on  sphere  diameter 

u  velocity  vector  of  the  fluid 

v  velocity  vector  of  the  sphere 

initial  maximum  induced  velocity  due  to 
the  vortex  tube  normalized  by  v0 
v0  initial  injection  velocity  of  the  sphere 

tf  dimensionless  time  normalized  by  a/v0 

t*  dimensionless  time  normalized  by  u 

X,  E  inertial  cylindrical  coordinates 

x,cr  noninert ial  cylindrical  coordinates 

X,  y,  Z  inertial  Cartesian  coordinates 

X'0 ,  Yc'  initial  particle  position 
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normalized  by  the  particle  radius 
noninertial  Cartesian  coordinates 

Greek  symbols 

pi 

fluid  density 

pp 

particle  density 

PS 

fluid  viscosity 

vj 

fluid  kinematic  viscosity 

nonorthogonal  generalized  coordinates 

(Pc 

initial  radius  of  vortex  core 

u' 

dimensionless  frequency  normalized  by  a/ 

Superscript 

/ 

nondimensional  quantity 

Subscript 

/ 

carrier  fluid 

0 

initial  quantity 

P 

spherical  particle 

1.  Introduction 


Accurate  prediction  of  particle  (or  droplet)  disper¬ 
sion  is  important  in  many  turbulent  two-phase  flows 
such  as  spray  combustion  and  atmospheric  dispersion 
of  pollutants.  Since  a  general  analytical  solution  of  the 
Navier-Stokes  equations  for  a  flow  around  a  sphere  is 
not  possible,  only  a  numerical  solution  of  these  equa¬ 
tions  can  provide  accurate  information  about  the  flow 
field.  The  forces  on  the  particle  can  then  be  computed 
by  integrating  the  normal  and  shear  stresses  around  the 
particle,  and  then  Newton’s  law  is  applied  to  obtain  the 
acceleration  of  the  particle.  However,  since  these  equa¬ 
tions  are  unsteady  and  three-dimensional,  they  require 
excessive  computing  time.  It  is  not  possible  to  use  this 
method  to  predict  the  simultaneous  motion  of  many 
particles  in  a  typical  two-phase  flow  with  the  present 
and  foreseeable  computing  capabilities. 

The  simplest  method  to  describe  the  forces  on  a  par¬ 
ticle  would  be  using  the  steady  drag  from  the  standard 
drag  curve  plus  the  net  gravitational  force  on  the  par¬ 
ticle  according  to 

mP~j:  —  icD,tdTra2pj  \  u  -  v  |  (u  -  »)  +  (mp-m/)g 

(1) 
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where  Costd  is  the  drag  coefficient  from  the  (steady) 
standard  drag  curve,  u  is  the  carrier  fluid  velocity,  v 
is  the  sphere  velocity,  a  is  the  sphere  radius,  and  m/ 
is  the  mass  of  the  fluid  displaced  by  the  sphere.  Many 
application-oriented  studies  of  dilute  particle  dispersion 
are  based  on  Equation  (1).  In  the  case  of  unsteady  flow, 
this  relationship  is  an  approximation  that  can  be  valid 
only  if  the  response  time  of  the  particle  is  much  larger 
than  that  of  the  flow. 

Several  equations  accounting  for  the  unsteadiness  of 
the  particle  motion  have  been  developed  wherein  a  su¬ 
perposition  of  the  steady  drag  and  the  unsteady  (his¬ 
tory)  drag  is  used  to  obtain  the  forces  on  the  particle. 
The  available  particle  equations  are  reviewed  next. 

Basset1,  Boussinesq2,  and  Oseen3  studied  the  un¬ 
steady  rectilinear  motion  of  a  sphere  in  a  stagnant  in¬ 
compressible,  viscous  fluid.  They  solved  the  Navier- 
Stokes  equations  for  a  creeping  flow  by  neglecting  the 
advective  acceleration  terms  and  derived  the  following 
expression  for  the  acceleration  of  the  sphere: 


dv  1  dv  .  . 

mp—  =  —6irafijv  -  - mj —  +  (mp  -  mj)g 


-  6a? y/WJiJpJ  (  ^j= L  dr  (2) 


Maxey  and  Riley4  re-derived  from  first  principles 
the  following  equation  for  the  motion  of  a  sphere  in 
a  nonuniform  creeping  flow: 


mP^  =  6*anf(u- 


x  1  d(u  —  v)  Du 

+  2”1'  it  +  m,~sr 


+  6aiy/WitJpJ  j  ^-j=JZ-dr  +  (mp-mrfg,  (3) 


where  the  Faxen  forces  are  not  shown  here. 

Based  on  the  above  two  expressions  (2)  and  (3),  one 
may  write  an  equation  for  the  motion  of  a  sphere  for 
finite  Reynolds  number: 


dv 

^Tt 


CD,td*a?Pj  |  u  -  v  |  (u  —  v)  + 

d{u  —  v)  Du 

dt 

*  d{u  —  v)/dr 


«mj-  +  m}—  +  (m„ 


-  mj)g  + 

.  2  _ _  [l  d(u  —  v)/dr  ... 

6a?  y/WJTjpj  — dr,  (4) 


where  the  first  term  is  an  empirical  modification. 

Odar  and  Hamilton5  and  Odar6  experimentally  ex¬ 
amined  the  force  on  a  guided  sphere  rectilinearly  oscil¬ 
lating  in  an  otherwise  stagnant  fluid  for  0  <  Re  <  62. 
They  proposed  an  equation  for  the  motion  of  a  sphere 
in  a  flow  of  finite  Reynolds  number  based  on  their  ex¬ 
perimental  study  as: 


dv 


ms 


dt 


-)^CDstdTra?pj  |  v  |  v  -  Ca ~ 
Ch  6a2y/¥JTfpJ  f  dr  ,  (5) 

J o  y/t  —  T 


with  Ca  and  Ch  obtained  experimentally  and  given  by 

Ca  =  2.1  -  O.mMh/il  +  O.UMh) 

Ch  =  0.48  +  0.52M^1/(1  +  Mai)3  , 


where  MA\  is  the  dimensionless  acceleration  defined  by 
2a 


MAl  = 


i« 


-«|2 


1^1- 


Mei  et  al7  studied  an  unsteady  flow  over  a  station - 
ary  sphere  with  small  fluctuations  in  the  freestream 
velocity  at  finite  Reynolds  number  (0.1  <  Re  <  40) 
using  a  finite  difference  method  and  found  that  the 
Basset-force  term  in  the  equation  of  particle  motion 
should  have  a  kernel  which  must  decay  much  faster 
than  l/y/t  —  t  at  large  time.  Mei  and  Adrian8  and 
Mei9  considered  the  same  problem  as  Mei  et  al.7  but 
fox  St  «  Re  «  l  using  a  matched  asymptotic  expan¬ 
sion,  where  St  is  Strouhal  number  based  on  the  angu¬ 
lar  frequency  of  the  freestream  and  the  sphere  radius. 
They  proposed  a  modified  expression  for  the  Basset- 
force  term  on  the  basis  of  the  analytical  result  at  small 
Reynolds  number  for  low  frequency,  the  numerical  re¬ 
sult  at  finite  Reynolds  number  for  low  frequency,  and 
the  unsteady  Stokes’  result  for  high  frequency.  Their 
proposed  equation  is: 

dv  1  _  o  i  .  /  x 

mp—  =  -CDstd™  pj\u~v\(u-v)  + 

1  (Du  dv 

2mj  \Dt  dt 

6irpja  f  K(i  —  T>T)~~i — —dr  + 

J -co 

(mp  -mf)g,  (6a) 
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+  mj 


Du 

Dt 


with  the  broad-frequency-range  approximation  for  the 
integral  kernel  given  by 


K(t  -  t,  t)  = 


{ 


r  u  \  0.25 

K(t  —  T)l/f 

a 2 


+ 


>|»(r)-r(T)|3 

.2 


where  fff(Rct)  =  0.75  +  0.105iZe*(r);  Ret  =  |tt(r)  — 
v(r)\2a/vj. 

Maxey10  included  the  effect  of  the  initial  velocity  dif- 
ference  between  the  sphere  and  the  carrier  fluid  in  the 
particle  motion  equation  of  Maxey  and  Riley.4  The  ad¬ 
ditional  term  is:  6a2y/irpjpf(u(0)  —  v(0))/y/t. 

In  Section  2,  the  numerical  solutions  from  the  above 
equations  will  be  compared  with  those  from  the  full 
Navier-Stokes  equations  for  unsteady,  axisymmetric 
flow  around  a  freely  moving  sphere  initially  injected 
into  an  oscillating  fluid,  and  a  new  equation  for  arbi¬ 
trary  rectilinear  particle  motion  is  proposed.  In  Section 
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3,  the  numerical  solutions  from  the  above  equations 
will  be  compared  with  those  from  the  full  Navier-Stokes 
equations  for  unsteady,  three-dimensional  flow  around 
a  freely  moving  sphere  interacting  with  a  large  vortex 
tube,  and  a  new  equation  for  arbitrary  two-  and  three- 
dimensional  particle  motion  will  be  proposed.  Section 
4  provides  a  summary  of  the  work.  The  gravity  force  is 
neglected  in  the  computation  in  the  following  sections. 


2.  Unsteady,  axisymmetric  flow  around  a  sphere 
injected  into  an  oscillating  fluid 

2.1  Flow  description 

Consider  an  unsteady,  axisymmetric,  incompressible, 
laminar  flow  generated  by  a  spherical  particle  injected 
into  a  constant  property  Newtonian  fluid  oscillating 
with  time  in  the  same  direction  as  that  of  the  parti¬ 
cle  motion  as  shown  in  Figure  1.  The  net  gravity  force 
acting  on  the  particle  is  neglected.  The  origin  of  a  non- 
rotating  noninertial  reference  frame  is  chosen  at  the 
center  of  the  particle. 

Three  coordinate  systems  are  used  in  our  formula¬ 
tion:  the  inertial  (fixed  in  space)  cylindrical  coordinates 
(X,  E),  the  noninertial  cylindrical  coordinates  (x,<r), 
and  the  generalized  coordinates  (£,»?)•  The  origin  of 
the  coordinates  (x,  <r)  coincides  with  the  sphere  center. 
The  coordinate  £  gives  the  radial  and  rj  gives  the  an¬ 
gular  direction  with  respect  to  the  sphere  and  are  used 
for  the  numerical  solution  of  the  Navier-Stokes  equa¬ 
tions.  The  generalized  coordinate  system  can  be  easily 
adapted  to  two-dimensional  or  axisymmetric  arbitrary 
geometries. 

The  base  flow  in  the  far  field  oscillates  with  time  in 
X-direction  and  is  expressed  as 

ux(t)  =  ai|t/|sinwf 
Uc(t)  =  0  , 

where  ai  is  a  constant  controlling  the  amplitude  and 
u  is  the  angular  frequency.  The  associated  far  pressure 
field  can  be  obtained  from  the  Navier-Stokes  equations 
and  is  given  as 

Pb(x,t)  =  — sinutf  +  \v\u>  cosu>t)X  +  pref  , 

where  prej  is  the  reference  pressure  at  X  =  0.  X  is 
related  to  x  as  X  =  x  +  Xp ,  where  Xp  is  the  distance 
travelled  by  the  particle  and  measured  from  the  origin 
of  the  inertial  coordinates  (X,  E). 

2.2  Governing  equations  and  boundary  condi¬ 
tions 


«(£  +  f +v-vf)=-v,  +  „v’v.  <») 

The  governing  equations  (7)  and  (8)  are  are  nondi- 
mensionalized  using  the  sphere  radius  a  as  the  charac¬ 
teristic  length  and  the  initial  injection  velocity  of  the 
sphere  vQ  as  the  characteristic  velocity.  The  nondimen- 
sionalized  equations  are  cast  in  terms  of  the  generalized 
coordinates  (f ,  rj)  to  treat  an  axisymmetric  body  of  ar¬ 
bitrary  shape.  The  numerical  integration  is  performed 
using  a  cubic  computational  mesh  with  equal  spacing 
(6£  =  6t)  =  1).  In  the  present  study,  a  spherical  do¬ 
main  is  used,  and  the  grid  reduces  to  an  orthogonal, 
spherical  grid.  The  grids  are  denser  near  the  surface  of 
the  spherical  particle,  and  the  grid  density  in  the  radial 
direction  is  controlled  by  the  stretching  function  devel¬ 
oped  by  Vinokur.11  The  domain  of  the  flow  is  bounded 
by  I  <  £  <  Nj,  1  <  rj  <  N2,  where  £  =  1  and  N\ 
correspond,  respectively,  to  the  sphere  surface  and  the 
farfield  boundary  surrounding  the  sphere;  rj  =  1  and  N2 
denote,  respectively,  the  positive  x-axis  (downstream) 
and  the  negative  x-axis  (upstream). 

The  velocities  on  the  sphere  surface  are  zero  due  to 
the  no-slip  condition,  and  the  pressure  on  the  sphere  is 
obtained  from  the  momentum  equation.  The  boundary 
conditions  are 


dp 

dn 


d2Vn  dvn 
**  dn 2  PS  dt  ’ 


Vi  =  K,  =  0at£  =  l  (9a) 


v  =  ph,  V-  -  n.  -  v.  Ve  =  0  at  (  ~  Ar;„.  ^  ”  <  N2 

(9b) 


P  ~  Pb, 


dK 

dx 


^  =  0atf  =  JV1,l<ii< 

Ox 


N2m 

(9c) 


dp 

drj 


— -  =  0,  Vc  =  0  at  r)  =  1  and  N2 
dr) 


(9d) 


where  Vx  and  Va  are  the  fluid  velocities  with  respect  to 
the  sphere  in  the  x  and  a  directions,  respectively,  Vn 
is  the  fluid  velocity  with  respect  to  the  sphere  in  the 
direction  normal  to  the  sphere  surface,  and  vn  is  the 
sphere  velocity  in  the  direction  normal  to  the  sphere 
surface,  n  denotes  the  direction  normal  to  the  sphere 
surface,  d/dn  =  \/£x  +  and  rj  =  N2m  denotes 

the  mid-plane  between  rj  —  1  and  N2- 

In  order  to  start  the  numerical  solution  of  equations 
(7)  and  (8),  we  provide  initial  velocity  by  superposing 
the  initial  velocities  of  the  base  flow  and  the  sphere, 
and  the  no-slip  condition  on  the  sphere  surface: 

Po  =  Pb ,  vx0  =  -Vo,  Vac  =  o,  except  at  (  =  1 

(10a) 


The  continuity  and  momentum  equations  to  be 
solved  are: 

V-V  =  0  (7) 
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Po  =  Pb ,  V„  =  Vao  =  0  at  £  =  1,  (10b) 


2.3  Numerical  solution 

The  three-dimensional  algorithm  to  be  used  in  pre¬ 
dicting  the  unsteady,  three-dimensional  flow  around 
a  sphere  interacting  with  a  large  vortex  (Section  3) 
has  been  described  in  References  12  and  13.  Here, 
an  axisymmetric,  implicit,  finite- difference  algorithm 
has  been  developed  to  solve  simultaneously  the  set 
of  the  discretized  partial  differential  equations.  The 
method  is  based  on  an  Alternating-Direction-Predictor- 
Corrector  (ADPC)  scheme  to  solve  the  time-dependent 
Navier-Stokes  equations.  ADPC  is  a  slight  variation  of 
Alternating-Direction-Implicit  (ADI)  method  and  im¬ 
plemented  easily  when  embedded  in  a  large  iteration 
scheme.  The  control  volume  formulation  is  used  to  de¬ 
velop  the  finite-difference  equations  from  the  governing 
equations  with  respect  to  the  generalized  coordinates 
(£,77).  An  important  part  of  solving  the  Navier-Stokes 
equations  in  primitive  variables  is  the  calculation  of  the 
pressure  field.  In  the  present  work,  a  pressure  correc¬ 
tion  equation  is  employed  to  satisfy  indirectly  the  con¬ 
tinuity  equation.  The  pressure  correction  equation  is  of 
the  Poisson  type  and  is  solved  by  the  Successive-Over- 
Relaxation  (SOR)  method. 

The  overall  solution  procedure  is  based  on  a  cyclic 
series  of  guess-and-correct  operations.  The  velocity 
components  are  first  calculated  from  the  momentum 
equations  using  the  ADPC  method,  where  the  pressure 
field  at  the  previous  time  step  is  employed.  This  esti¬ 
mate  improves  as  the  overall  iteration  continues.  The 
pressure  correction  is  calculated  from  the  pressure  cor¬ 
rection  equation  using  the  SOR  method,  and  new  es¬ 
timates  for  pressure  and  velocities  are  obtained.  This 
process  continues  until  the  solution  converges  at  each 
time  step. 

We  now  test  the  accuracy  of  the  solution  procedure 
by  predicting  the  axisymmetric  flow  over  a  solid  sphere. 

Here  we  examine  the  flow  generated  by  an  impul¬ 
sively  started  solid  sphere  in  a  quiescent  fluid  at  two 
Reynolds  numbers:  20  and  100.  The  time-dependent 
solution  converges  asymptotically  to  a  steady-state 
which  is  in  good  agreement  with  the  available  experi¬ 
mental  data  and  correlations  as  shown  in  Table  1.  Table 
1  lists  the  drag  coefficient  as  a  function  of  the  compu¬ 
tational  grid  density  at  Reynolds  numbers  20  and  100 
respectively,  and  compares  them  with  the  correlations 
of  Clift  ei  a/.14  Table  1  also  shows  the  separation  angle 
measured  from  the  front  stagnation  point  as  a  func¬ 
tion  of  grid  density  at  Reynolds  number  20  and  100, 
in  comparison  with  the  data  of  Taneda15  and  also  with 
the  correlations  of  Clift  et  a/.14  The  calculations  were 
performed  for  four  different  grids,  (N1XN2)  —  (21x21), 
(31  x  31),  (41  x  41),  and  (51  x  51)  in  a  computational 
domain  with  an  outer  boundary  located  at  21  sphere 
radii  from  the  sphere  center.  The  51  x  51  grid  is  used 
in  the  following  calculations. 

We  tested  the  solution  procedure  by  varying  the  far- 


field  boundary  condition  and  by  changing  the  location 
of  the  outer  boundary.  In  the  first  test,  the  far-field  out¬ 
flow  boundary  condition  was  changed  from  d<j)/dx  =  0 
{<j>  ~  Vx  and  Va).  to  d<f>/dr  —  0  There  was  almost  no 
difference  in  the  drag  coefficient  and  the  near  wake  size 
(the  separation  angle  and  length  of  the  recirculation 
eddy)  at  Reynolds  numbers  20  and  100.  Our  calcula¬ 
tion  shows  that  separation  does  not  occur  at  Reynolds 
number  20.  In  the  second  test,  the  location  of  the  outer 
boundary  in  downstream  was  changed  from  21  to  41 
sphere  radii.  There  was  virtually  no  change  in  the  drag 
coefficient  and  the  near  wake  size  at  both  Reynolds 
numbers. 

2.4  Comparison  of  the  solution  of  the  previous 
equations  with  that  of  the  Navier-Stokes 
equations 

We  now  compare  the  numerical  solutions  from  the 
equations  introduced  in  Section  1  with  those  from  the 
full  Navier-Stokes  equations  for  unsteady,  axisymmetric 
flow  around  a  freely  moving  sphere  initially  injected  into 
a  still  fluid  (the  case  of  w  =  0  in  Section  2.1). 

Figure  2  shows  the  drag  coefficients  of  the  sphere  as 
a  function  of  time  (0  <  t*  <  200)  with  initial  particle 
Reynolds  number  Rei0  =  150  and  the  sphere/fluid  den¬ 
sity  ratio  pr  =  5.  At  t'  =  200,  the  particle  Reynolds 
number  reduces  to  4.  The  Basset  history  term  in  Equa¬ 
tions  (4)  and  (5)  causes  too  low  a  value  for  the  drag 
coefficient  compared  with  the  Navier-Stokes  solution. 
The  drag  coefficient  from  Equation  (6a)  proposed  by 
Mei  and  Adrian  (1992)  is  the  closest  to  that  from  the 
Navier-Stokes  equations,  but  with  increasing  discrep¬ 
ancy  as  if  increases. 

Figure  3  shows  the  drag  coefficients  of  the  sphere  as 
a  function  of  time  (0  <  t*  <  400)  with  the  same  particle 
Reynolds  number  as  in  Figure  2  but  with  the  density 
ratio  pr  =  200.  At  t!  =  400,  the  particle  Reynolds  num¬ 
ber  becomes  84.  The  Basset  history  term  in  Equations 
(4)  and  (5)  still  results  in  lower  values  for  the  drag  coef¬ 
ficient.  Now,  the  solution  of  Equation  (6a)  gives  a  good 
approximation  to  that  of  the  Navier-Stokes  equations. 
However,  it  is  noted  that  the  solution  of  Equation  (1) 
also  gives  a  good  approximation  and  is  very  close  to 
that  of  the  Navier-Stokes  equations.  Thus,  it  is  seen 
that  for  high  density  ratio  ( pT )  the  deviation  from  the 
Navier-Stokes  solution  is  reduced. 

2.5  New  improved  equation  for  particle  motion 

We  now  investigate  the  restrictions  imposed  on  the 
derivation  of  Equation  (6a)  proposed  by  Mei  and 
Adrian8  and  propose  an  improved  equation  for  the  par¬ 
ticle  motion. 

Two  restrictions  were  imposed  on  the  derivation  of 
Equation  (6a).  First,  the  integral  kernel  was  developed 
under  the  assumption  of  small  amplitude  oscillation  of 
the  free  stream.  Secondly,  the  effect  of  the  drag  due 
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to  the  initial  relative  velocity  between  the  particle  and 
the  carrier  fluid  was  neglected. 

As  seen  in  the  previous  section,  Equation  (6a)  pro¬ 
duces  a  better  solution  for  the  drag  of  a  sphere  with 
higher  density  ratio  than  for  the  case  with  low  density 
ratio;  in  other  words,  Equation  (6a)  produces  a  better 
solution  for  the  drag  of  a  sphere  with  lower  deceleration 
than  with  higher  deceleration. 

The  above  analysis  indicates  the  need  for  the  sphere 
motion  equation  to  account  for  the  initial  velocity  dif¬ 
ference  and  to  have  an  integral  kernel  weighted  by  the 
acceleration  magnitude.  The  weighting  function  con¬ 
tains  the  time  derivative  of  the  relative  velocity  Mai 
and  the  ratio  <j>r  of  Ma2  to  Mai-  Mai  was  defined 
before  in  Equation  (5)  and  Mai  and  <j>r  are  defined  by 


where  Ret  =  |u(r)  —  v{r)\2a/vj  and  Ret0  =  |u(0)  — 
v(0)\2a/i/j . 

The  six  constants  c,*  ( i  =  1,  ..,6)  in  the  above  equa¬ 
tions  are  determined  by  comparing  the  solutions  from 
Equation  (11a)  with  those  from  the  Navier-Stokes  equa¬ 
tions  and  given  for  0  <  u/'  <  oo  by 

ci  =  2.5,  C5  =  0.126,  ce  =  15, 

c2  =  45.5,  c3  =  0.03,  c4  =  0.1 .  (llh) 

For  low  frequencies  a/  <  0.3,  the  following  values  of 
c2,  c3,  and  c4  provide  slightly  better  results  than  those 
of  Equation  (llh).  The  values  of  ci,  C5,  and  c$  are  kept 
fixed  as  in  (llh). 

c2  =  13.9,  c3  =  0.12,  c4  =  0.5  .  (lli) 


MA2(t)  = 


(2a)2 

|«  - t>|3 


<p\u  —  v\ 

dp 


I ;  Mf) 


Ma2 

Mai 


These  dimensionless  groups  can  be  introduced  through 
dimensional  analysis  to  obtain  the  forces  on  the  particle 
of  unsteady  motion. 

Now,  we  propose  the  following  equation  for  particle 
motion. 

Tn?~Jl  =  \cDstd*a2pf  |  «  -  v  |  (u  -  t>)  + 


6717*/ a  Ki(t)  [u(0)  -  v(0)]  + 

(mp  -  mj)g  (11a) 

with  the  integral  kernel  K(t  —  r,  r)  given  by 


K(t-r,  r) 

G(r) 

G(r) 

P 

!h 


■II 

V(f  —  r)vj ’ 

a3  j 

0.5/ Cj 

+ 

n  |u(r)  —  t)(r)|3  . 
.2  avjfjj(Ret)  K 

~rf 

1/ci 

1 

1 

+  /?  s/Mai1 

F) 

C2 

T 

+  rfr#7[es(rfr  +  #4)] 

=  0.75  +  c$Ret(r) , 


(lib) 

(11c) 

(lid) 

(lie) 


When  G(r)  equals  unity,  the  present  integral  kernel 
(Equation  (lib))  is  similar  to  Equation  (6b)  by  Mei 
and  Adrian,8  but  the  values  of  c\  and  C5  in  the  present 
integral  kernel  are  different  from  those  used  in  Equa¬ 
tion  (6b).  The  values  2  and  0.105  in  Equation  (6b)  cor¬ 
responding  to  c\  and  C5  were  determined  respectively  by 
an  interpolation  and  a  curve  fitting  (Mei  and  Adrian8). 

Equation  (11c)  shows  that  as  Mai  is  reduced, 
G(r)  approaches  unity,  and  the  present  kernel  (Equa¬ 
tion  (lib))  becomes  similar  in  form  to  the  kernel  of 
Equation  (6b).  On  the  other  hand,  as  Mai  becomes 
large,  G(r)  approaches  zero,  and  the  present  kernel 
becomes  of  the  same  form  as  that  of  Basset1  (Equa¬ 
tion  (4)).  The  (3  in  the  expression  of  G(r)  is  not  a  con¬ 
stant  but  a  function  of  <f>r,  the  ratio  of  Ma2  to  Mai- 
This  function  behaves  as  follows.  (3  ~  c2(  1  —  <£r/c3)  for 
<j>r  «  1,  and  /?  ~  c2c3/<^4  for  <j>r  »  1.  Also,  it  can 
be  shown  that  G(r*)  ~  u/"°’4  as  a/  »  1  when  r  is 
normalized  by  u  (r*  =  ru). 

G\  (Equation  (llg))  is  obtained  from  G(r  = 
0)  (Equation  (11c))  as  follows.  From  the  equa¬ 
tion  of  particle  motion  Equation  (11a),  it  can  be 
shown  that  dv/dt  ~  fcf”0,5  as  t  — ►  0,  where  k  = 
(4.5/a)  y/vj  fic  (pr  +  O.S^MO)  -  t>(0)).  Employing 
this  form  of  dv/dt ,  the  expressions  of  Mai,  Ma2,  and 
<j>r  as  t  — ►  0  can  be  also  derived.  For  example,  <t>r  is  ex¬ 
pressed  by  <j>r  ~  ar’1/(u(0)  -  u(0)).  Finally,  G(r  — ►  0) 
is  obtained  as 

G(T  0)  “  1  +  0^4-0.25^-0.85^  +  0.5)-0  5 


and  with  the  function  K\{i)  given  by 


Ki(t) 


Cx 


{[^P 


+ 


p  K0)~t>(0)|3  211/Cl) 
[2  avjpH{Ret0)  J  J 


1  +  ceRefo  2S(Pr  +  0.5)-0-5  ’ 


(Ilf) 

(Hg) 


where  a  =  c2c3(9/>/7r)o'5(|ti(0)  —  v(0)|/a)C4“°‘25.  This 
equation  shows  that  G(0)  =  0  when  c4  <  0.25  but 
G(0)  =  1  when  c4  >  0.25.  However,  the  drag  coefficient 
from  the  numerical  computation  of  Equation  (11a)  with 
Gi  =  G(0)  =  0  or  Gi  =  G(0)  =  1  is  too  high  or  low 
for  some  initial  period  compared  with  that  from  the 
Navier-Stokes  equations.  The  drag  coefficient  from  the 
numerical  computation  of  Equation  (11a)  with  G 1  = 
G(r=0,  c4=0.25)  is  also  too  low  for  some  initial  period 
compared  with  that  from  the  Navier-Stokes  equations. 
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Therefore,  G\  was  introduced  as  Equation  (llg)  by  fol¬ 
lowing  the  form  of  G(r  — ►  0),  and  the  coefficient  ce  in 
the  expression  of  G\  is  determined  by  numerical  opti¬ 
mization. 

Figure  4  shows  the  drag  coefficients  of  the  sphere  as  a 
function  of  time  for  the  same  conditions  as  in  Figure  2, 
but  here  the  drag  coefficients  were  computed  from  the 
new  proposed  equation  (11a)  and  the  previous  equa¬ 
tions  including  an  initial  velocity  difference  term.  We 
note  that  the  initial  velocity  difference  terms  are  differ¬ 
ent  among  the  equations  because  the  integral  kernels 
are  different  among  the  equations.  The  initial  veloc¬ 
ity  difference  term  added  to  the  BBO  equation  (Equa¬ 
tion  (4))  is  the  term  given  by  Maxey10  and  shown  in 
Section  1.  The  initial  velocity  difference  term  added  to 
Equation  (6a)  by  Mei  and  Adrian8  is  the  term  given  by 
setting  G\  =  1  with  c\  =  2  and  C5  =  0.105  in  Equa¬ 
tions  (Ilf)  and  (lie).  It  is  shown  in  Figure  4  that  the 
present  equation  of  sphere  motion  gives  the  best  solu¬ 
tion. 

Comparing  Figures  2  and  4,  we  see  that  the  initial 
velocity  difference  term  with  f-0,5  improves  Equation 
(4).  However,  the  appropriate  decay  of  this  term  is  t~ 2 
for  finite  initial  particle  Reynolds  number  (see  Equa¬ 
tion  (Ilf)).  However,  it  is  interesting  to  note  that  the 
initial  velocity  difference  term  decays  as  5  when 
Ret0  — ►  0  (see  Equations  (llg)  and  (Ilf)). 

Figure  5  shows  the  drag  coefficients  of  the  sphere  as 
a  function  of  time  (0  <  V  <  50)  for  the  same  conditions 
as  in  Figure  4  except  that  Reto  =  38.  At  V  =  50,  the 
particle  Reynolds  number  becomes  2.18.  It  is  shown 
that  the  new  equation  (11a)  of  sphere  motion  gives  the 
best  solution  for  low  initial  particle  Reynolds  number 
as  well. 

Figures  6  and  7  show  the  drag  coefficients  of  the 
sphere  as  a  function  of  time  (0  <  V  <  200)  with 
Reto  =  150  and  the  density  ratio  pr  =  5.  The  base  flows 
in  Figures  6  and  7  oscillate  with  a/  =  0.4  and  ot\  —  0.02, 
and  u /  =  0.8  and  ai  =  0.01,  respectively.  Both  figures 
show  that  Equation  (6a)  with  an  initial  velocity  differ¬ 
ence  term  produces  higher  drag  coefficients  (except  for 
some  initial  period)  than  do  the  Navier-Stokes  equa¬ 
tions.  Again,  the  new  equation  (11a)  produces  very 
good  agreement  with  the  Navier-Stokes  equations. 

Figures  8  shows  the  drag  coefficients  of  the  sphere 
as  a  function  of  time  (0  <  tf  <  200)  with  Ret0  =  150 
and  the  density  ratio  pr  =  200.  The  base  flow  oscil¬ 
lates  with  u/  =  0.1  and  cti  =  0.18.  This  figure  shows 
that  Equation  (11a)  produces  slightly  better  drag  co¬ 
efficient  than  does  Equation  (6a)  with  an  initial  veloc¬ 
ity  difference  term  compared  with  the  drag  coefficient 
from  the  Navier-Stokes  equations.  The  better  perfor¬ 
mance  of  Equation  (6a)  for  the  case  of  higher  density 
ratio  is  due  to  small  dimensionless  acceleration  Max  in 
the  case  of  higher  density  ratio  as  shown  in  Figure  9, 
which  shows  Mai  as  a  function  of  time  for  the  cases  of 
pr  —  h  (with  u/  ~  0.1  and  c*i  =  0.06)  and  pr  =  200 


(Figure  8).  When  Max  becomes  small,  the  function 
G(t)  in  Equation  (lib)  approaches  unity,  and  thus  the 
integral  kernel  of  Equation  (11a)  approaches  that  of 
Equation  (6a). 

Figure  10  shows  the  drag  coefficients  of  the  sphere  as 
a  function  of  time  (0  <t*<  200)  for  the  same  condi¬ 
tions  as  in  Figure  8  except  that  the  drag  coefficients  are 
obtained  from  different  equations.  Neglecting  the  his¬ 
tory  term  and  the  other  terms  in  Equation  (11a)  causes 
a  phase-lag  to  the  drag  coefficient. 

3.  Unsteady,  three-dimensional  flow  around 
a  sphere  interacting  with  a  large  vortex 

3.1  Flow  description 

We  consider  the  time-dependent,  three-dimensioned, 
incompressible,  flow  around  a  small  spherical  solid  par¬ 
ticle  injected  into  a  counter-clockwise  rotating  large 
vortex  which  is  located  at  the  origin  of  the  coordinates 
pf,y,  Z)  fixed  in  space  as  shown  in  Figure  11.  The 
net  gravity  force  acting  on  the  particle  is  neglected. 
The  origin  of  a  nonrotating  noninertial  reference  frame 
(x,y,  z)  is  chosen  at  the  center  of  the  particle. 

The  initial  velocity  field  induced  by  the  vortex  tube 
is  analytically  computed  by  considering  the  evolution 
of  a  point  vortex  and  is  given  by 

-“'(-4^rk)>i  <m> 

,  .  d 

0  Avs  (L12)2  ’ 

where  R  =  y/X2  -f  Y2,  <j>  =  arctanY/X,  and  t0  is  a 
parameter  defined  by  initial  size  of  the  vortex  core  (<r0) 
and  the  fluid  kinematic  viscosity 

The  initial  maximum  induced  velocity  umoxo  due  to 
the  vortex  tube  occurs  on  the  edge  of  the  vortex  core 
R  =z  <r0  at  t  =  0  and  can  be  obtained  from  Equar 
tion  (12a).  The  pressure  field  due  to  the  vortex  tube 
is  also  analytically  computed  by  integrating  the  ra¬ 
dial  component  of  the  momentum  equation  which  is 
pUl/R  =  dp/dR. 

It  is  assumed  that  the  vortex  is  so  large  that  the 
flow  field  induced  by  the  vortex  is  not  affected  by  the 
moving  sphere  except  locally  in  the  region  around  the 
sphere.  The  boundary  conditions  of  the  flow  field  are 
obtained  by  superimposing  the  sphere  velocity  and  the 
induced  velocity  due  to  the  vortex  tube  at  the  com¬ 
putational  outer  boundary  relative  to  the  sphere,  and 
the  pressure  field  due  to  the  vortex  tube  is  also  im¬ 
posed  at  the  computational  outer  boundary.  Although 
the  flow  computation  is  three-dimensional,  the  particle 
path  remains  in  the  plane  of  symmetry  (X  -  Y  plane). 
Therefore,  the  trajectory  computation  is  made  for  the 
case  of  a  freely  moving  particle  in  two-dimensions.  Af¬ 
ter  computing  the  forces  on  the  sphere,  the  deceleration 
(or  acceleration)  of  the  sphere  is  obtained  via  Newton’s 
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second  law  of  motion,  and  then  the  new  location  of  the 
sphere  is  obtained. 

The  rotation  of  the  sphere  due  to  the  torque  on  the 
sphere  is  neglected  because  the  torque  is  very  small  as 
will  be  shown  in  Section  3.2. 

3.2  Equation  for  particle  motion  in  vector  form 

Figure  12  shows  the  lift  and  moment  coefficients  of 
the  sphere  as  a  function  of  time  (0  <  t*  <  600)  com¬ 
puted  by  solving  the  Navier-Stokes  equations  for  ini¬ 
tial  particle  Reynolds  number  Reto  =  108.2  and  par¬ 
ticle  position  =  (250,0)  with  density  ratio 

pr  =  pp/pj  —  200.  The  vortex  tube  has  initial  core 
radius  (a'0 )  equal  to  200  and  the  maximum  induced 
velocity  due  to  the  vortex  (n!maxo)  equal  to  2.  It  is 
found  that  the  lift  and  moment  coefficients  are  small 
( \Cl\  <  0-16  and  \Cm\  <  0.03)  compared  with  the  drag 
coefficient,  which  is  shown  in  Figure  13.  The  lift  and 
moment  coefficients  are  small  because  the  strain  rate 
of  the  flow  across  the  sphere  is  small  when  the  sphere 
size  is  much  smaller  than  the  vortex  tube  as  in  the  flow 
considered  here. 

Since  the  lift  and  moment  coefficients  are  small,  we 
postulate  that  the  force  due  to  the  history  term  is 
aligned  with  the  direction  of  the  relative  velocity  vector 
and  propose  the  following  equation  for  particle  motion 
in  vector  form. 


as  defined  in  Equations  (11c),  (lie),  and  (llg),  respec¬ 
tively.  For  rectilinear  motion,  Equation  (13a)  reduces 
to  Equation  (11a). 

The  six  constants  c,  (i  =  1, 2, 6)  in  the  above  equa¬ 
tions  are  the  same  as  those  in  Equations  (1  lh)  and 

(in). 

It  should  be  noted  that  Equation  (13a)  does  not  in¬ 
clude  terms  associated  with  the  lift  force  in  cases  of  a 
rotating  particle  or  a  particle  in  a  strong  shear  flow. 

Figure  13  shows  the  drag  coefficients  of  the  sphere 
as  a  function  of  time  computed  by  solving  the  Navier- 
Stokes  equations,  Equation  (13a),  and  Equation  (1) 
for  the  same  parameters  as  used  in  Figure  12.  It  is 
shown  that  the  drag  coefficient  from  Equation  (1)  is 
in  phase-lag  compared  with  that  of  the  Navier-Stokes 
equations.  Neglecting  the  history  term  and  the  other 
terms  in  Equation  (1)  causes  a  phase-lag  on  the  drag 
coefficient. 

Figure  14  shows  the  trajectory  of  the  sphere  in  the 
X  —  Y  symmetry  plane.  The  semi-circle  in  the  figure  is 
the  initial  vortex  core.  Equation  (1)  predicts  higher  Y- 
location  of  the  sphere  than  do  the  Navier-Stokes  equa¬ 
tions.  However,  the  deviation  of  the  trajectory  from 
Equation  (1)  is  small  compared  with  that  from  the 
Navier-Stokes  equations  because  the  dimensionless  ac¬ 
celeration  of  the  sphere  is  small  (Mai  <  0.0195)  due 
to  the  high  density  ratio  and  the  low  frequency  of  the 
large  vortex  tube. 


m. 


^■7  =  \cD,tdTta2pj  |  u  -  v  |  (u  -  v)  + 


dt 


Du 

~Dt 


dv\  Du 

--di)+m>-n  + 

(u  —  v)  f*  ,  d|ti  —  v\ 

l  K{t-T’T)—sr- 

67Tfija  Ki(t)  [u(0)  —  v(0)]  + 

(mp  -  mf)g 

with  the  integral  kernel  K(t  —  r,  r)  given  by 


dr  + 


(13a) 


4.  Conclusion 

An  improved  equation  for  rectilinear  particle  motion 
has  been  proposed  which  includes  a  modified  history 
term  and  a  drag  force  due  to  the  initial  velocity  differ¬ 
ence  between  the  particle  and  the  carrier  flow. 

An  extension  of  the  equation  to  a  vector  form  for 
two-  or  three-dimensional  motion  has  been  also  pro¬ 
posed;  the  first  comparison  for  particle  motion  in  two 
dimensions  is  favorable.  More  computations  are  under¬ 
way  to  examine  the  equation  in  detail. 


K(t-r,  t) 
G(r) 


/  pr (t-r)v/ 
=  {[  a> 


1  0.5/ Ci 


x  |u(r)  -  t>(r)|3  _  2] 1/Cl 

1.2  aVjfl{Ret)  K  ’ 


—  Ci 

►  (13b) 


and  with  the  function  K\(t)  given  hy 


M"  * 

p  |tt(0)  -p(0)|3  21 
1  l2  ai/jffj(Reto)  \ 


where  Ret  =  |m(t)  —  v(t)\2ol/vj\  Ret0  —  lM(0)  — 
v(0)\2a/vj  and  G{t),  fjj,  and  G\  are  in  the  same  form 
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Ni  x  N2  x  jV3 

Cop 

Cov 

Cd 

C*D 

8. 

8\ 

Re  =  20 

21  x  21 

1.087 

1.789 

2.876 

166.3 

31  x  31 

1.057 

1.759 

2.816 

172.6 

41  x  41 

1.042 

1.733 

2.775 

180 

51  x  51 

1.038 

1.725 

2.763 

2.74 

180 

180 

Re  =  100 

21  x  21 

0.558 

0.590 

1.148 

124.1 

31  x  31 

0.533 

0.581 

1.114 

125.6 

41  x  41 

0.524 

0.580 

1.104 

126.2 

51  x51 

0.521 

0.580 

1.101 

1.09 

126.4 

126.5 

Table  1.  Drag  coefficient  and  separation  angle  as  a  function  of 
grid  density  at  Re  =  20  and  100,  where  *  denotes 
the  data  from  the  correlation  of  Clift  et  ai  [12] 
and  Taneda  [13]. 
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Figure  1.  Flow  geometry  and  coordinates 


Figure  2.  Drag  coefficients  as  a  function  of  time  obtained 
from  various  equations  for  Reio  =  150  with  pT  =  5. 
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Figure  8.  Drag  coefficients  as  a  function  of  time 
for  Re  to  =  150  and  u/  =  0.1  with  pr  —  200. 


Figure  10.  Drag  coefficients  as  a  function  of  time  for  the 
parameters  as  in  Fig.  8  but  with  different  equations. 
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The  Influence  of  an  Advecting  Vortex  on  the  Heat 
Transfer  to  a  Liquid  Droplet  * 

M.  Masoudit  and  W.  A.  Sirignano* 

Department  of  Mechanical  and  Aerospace  Engineering 
University  of  California,  Irvine 


Abstract 

The  three-dimensional  interaction  of  an  initially 
cylindrical  vortex  tube  with  a  droplet  in  a  uniform 
stream  is  investigated  through  a  numerical  solution 
of  the  Navier-Stokes  equations.  Particular  attention 
is  given  to  the  effect  of  the  vortex  on  the  droplet  heat 
transfer.  The  transient  response  of  the  droplet  Nus- 
selt  number  is  sensitive  to  the  geometrical  factors 
that  specify  the  vortex  initial  position  and  struc¬ 
ture.  The  time- averaged  Nusselt  number  is  approxi¬ 
mately  the  value  for  the  axisymmetric  case  when  the 
vortex  center  approaches  the  droplet  along  the  base 
flow  symmetry  axis.  Correlation  for  time- averaged 
values  of  Nusselt  number  is  reported.  The  widely- 
used  correlation  for  the  droplet  heat  transfer  in  an 
axisymmetric  flow  is  corrected  to  account  for  the  in¬ 
fluence  of  the  vortex.  Based  on  our  findings,  it  is 
speculated  that,  in  a  spray  combustion  system,  the 
vortex-droplet  interactions  within  the  Kolmogorov 
scale  can  have  significant  effects  on  the  droplet  con¬ 
vective  heat  transfer. 

Nomenclature 

a'  dimensional  droplet  radius, 

(characteristic  length) 

d  vortex  offset  distance  from  the  base  flow 

symmetry  axis  (normalized  by  a') 


Ni,N2,N3 

number  of  grid  points  in  (£,77,0 

Nu 

Nusselt  Number 

Pc 

Peclet  number 

Pr 

Prandtl  number 

P 

pressure 

// 

9 

heat  flux 

r,e,<t> 

spherical  coordinates 

Re 

base  flow  Reynolds  number  (based  on 

droplet  diameter) 

t  time  (normalized  by  a' /U^) 

T  temperature 
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u,v,w  flow  velocities  in  (x,  ?/,  z)  directions 

(normalized  by  U^) 
dimensional  free  stream  velocity, 
(characteristic  velocity) 

Vmax  maximum  tangential  velocity  of 

vortex  tube  (normalized  by  U^) 

V  velocity  vector 

x,  y,  z  Cartesian  coordinates 

Greek  symbols 

(£,  77,  C)  computational  coordinates 

cr  radius  of  vortex  tube  (normalized  by  af) 

T  vortex  tube  circulation 

ip  stream  function 

v*  kinematic  viscosity  of  the  gas-phase 

r  shear  stress 

Superscript 

'  dimensional  quantity 

Subscript 

o  initial  quantity 

ax  quantity  in  the  corresponding 

axisymmetric  flow  (no  vortex) 
g  quantity  in  the  gas-phase 

/  quantity  in  the  liquid  phase 

s  droplet  surface  (gas-liquid  interface) 

v  vortex  quantity 

1.  Introduction 

The  fluid  dynamics  and  heat  transport  for  a  cold  liq¬ 
uid  droplet  in  a  hot  gaseous  axisymmetric  environ¬ 
ment  is  a  well-understood  phenomenon  and  there 
exists  substantial  literature  exploring  many  differ¬ 
ent  aspects  of  such  problems1  .  There  however  is 
a  shortage  of  literature  exploring  droplet  heating 
and  vaporization  when  the  far-field  flow  embracing 
the  droplet  undergoes  temporal  and/or  spatial  vari¬ 
ations.  Existing  literature  has  focused  on  variations 
due  to  acoustical  waves2-4  .  Vortical  disturbances 
have  not  been  widely  examined. 
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Such  a  class  of  problems  appears  when  the  droplet 
transport  properties  are  subject  to  velocity  and  tem¬ 
perature  fluctuations  in  a  turbulent  flow,  such  as 
might  occur  in  a  liquid- fueled  combustor.  In  par¬ 
ticular,  in  a  spray-droplet  system,  the  droplet  size 
is  ~  100  microns.  In  the  turbulence  spectrum  for 
many  continuous  combustors,  this  length-scale  cor¬ 
responds  to  that  of  the  Kolmogorov  scale;  thus,  the 
droplet  transport  phenomena  can  be  subject  to  tur¬ 
bulent  effects  primarily  associated  with  those  of  the 
Kolmogorov  scale.  Moreover,  since  turbulence  could 
be  represented  as  a  manifestation  of  vortex  dynam¬ 
ics5  ,  it  is  useful  to  study  the  effect  of  an  array  of 
vortices  on  the  droplet  where  the  size  of  the  vortices 
is  comparable  with  that  of  the  droplet.  In  this  work, 
we  have  studied  the  effect  of  one  advecting  vortex  on 
a  droplet. 

2.  Flow  Description,  Governing  Equations, 

and  Vortex  Characteristics 

The  solution  for  the  velocity  field  in  this  problem  has 
been  reported  in  previous  publications6,7  .  Below, 
we  present  a  summary  of  our  approach  including  the 
governing  equations,  the  boundary  and  initial  condi¬ 
tions,  the  computational  approach,  and  the  reasons 
necessitating  a  parameter  study. 

Consider  the  three-dimensional,  unsteady  flow 
field  surrounding  and  within  a  cold  droplet  impul¬ 
sively  injected  in  a  hot  gaseous  environment  with  the 
droplet  subsequently  subjected  to  an  unsteady  inter¬ 
action  with  an  advecting  vortex  tube.  The  problem 
is  non-linear  by  nature.  Fundamental  fluid  dynamic 
aspects  such  as  lift,  drag,  and  moment  coefficients  of 
the  vortex-droplet  interaction  were  reported6,7  and 
thus,  here,  we  concentrate  on  our  recent  findings  on 
variations  in  the  droplet  heat  transfer.  Constant 
properties  are  assumed  in  both  the  gas  and  the  liquid 
domain  and  the  droplet  experiences  no  vaporization. 
The  deceleration  of  the  droplet  due  to  the  drag  force 
is  not  considered. 

Since  our  goal  is  to  study  the  flow  interaction  with 
a  liquid  droplet,  we  solve  the  governing  equations 
for  both  the  gas  and  the  liquid  phase.  These  are 
the  Navier-Stokes  and  the  thermal  energy  equations 
in  both  phases;  the  continuity  equation  is  satisfied 
through  pressure  correction.  The  equations  and  the 
boundary  conditions  are  non-dimensionalized  using 
the  droplet  radius  a*  as  the  characteristic  length,  the 
undisturbed  free  stream  velocity  as  the  charac¬ 
teristic  velocity,  and  the  ambient  gas  temperature 
Tg  as  the  characteristic  temperature.  The  governing 
equations  are: 


Gas  phase 


=  0 

DVg  „  2 

Dt  ~  VPg  +  Reg  V  V>' 


Dt  RCgPVg 


Liquid  phase 


v-v;  =  o 


Tsr  =  +  -kvlVi 


(1) 

(2) 

(3) 

(4) 

(5) 


DT,  2  „2m 
—  V  T) 

Dt  ReiPri  1 


(6) 


These  governing  equations  are  transformed  to  the 
coordinates  (f,77,C);  see  Fig-  1.  £  is  the  radial,  rj  is 
the  angular,  and  (  is  the  azimuthal  direction.  The 
numerical  integration  of  the  equations  is  performed 
using  a  computational  cubic  mesh  with  equal  spac¬ 
ing  (<5£  zz6tj  =  6(=:  1). 


Gas /Liquid  Interface  Conditions 

The  conditions  at  the  interface  are  based  on  the  prin¬ 
ciple  of  continuity  of  shear  stresses  (the  discontinu¬ 
ity  in  shear  stress  across  the  surface  due  to  surface 
tension  gradient  has  been  shown  to  be  negligible  in 
its  impact  on  droplet  heating),  zero  normal  veloc¬ 
ity,  continuity  of  tangential  velocities,  continuity  of 
the  heat  flux,  and  continuity  of  temperature.  Since 
the  interface  in  our  flow  is  always  spherical  (un¬ 
der  the  assumption  of  small  Weber  number),  these 
conditions  are  conveniently  cast  in  terms  of  spher¬ 
ical  coordinates  (r,  6 ,  <j>)  with  its  origin  at  the  cen¬ 
ter  of  the  droplet.  The  (£ ,  77,  Q  coordinates  have  the 
same  orientation  as  the  spherical  coordinates  (r,  6 ,  <j>) 
but  obey  an  imposed  stretching  allowing  a  relatively 
denser  grid  concentration  near  the  gas-liquid  inter¬ 


face  (droplet  surface). 

Tlyr0,s 

— 

Tg,r$,s 

Tl  ,r0,$ 

= 

T 9,r<f>t3 

Vt,e,. 

— 

Vi,*., 

= 

t,,s 

= 

Tg,s 

<l"s 

ZZ 

C 

Here,  Trg:,  and  rr4,iS 

are 

respectively  the  shear 

stresses  on  a  positive  r-plane  in  the  positive  6  and  <f> 
direction  and  qn  is  the  heat  flux  from  the  hot  am¬ 
bient  gas  into  the  cold  liquid  droplet.  The  interface 
condition  for  pressure  is  obtained  from  the  momen¬ 
tum  equation. 
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Gas-Phase  Boundary  Conditions 

(Ni,N2,Ns)  and  (Nu,  i\T2j  N$)  are  the  number  of 
grid  points  in  the  gas  and  liquid  domain,  respec¬ 
tively,  in  (£,77,C)  coordinates,  f  at  Nu  and  N\  are 
the  droplet  surface  and  the  gas  far-field,  respectively. 
The  imposed  far-field  pressure,  gas  velocities  in  the 
x,  y,  z  directions  and  gas  temperature  are 

p  —  0,  u  =  v  =  0,  w  —  T  —  \  at  f  =  TVi  and 
N2mid  2  (upstream) 


Next,  we  use  the  interface  conditions  to  solve  for 
the  liquid-phase  boundary  values,  followed  by  the  se¬ 
quential,  iterative  solution  of  the  liquid-phase  equa¬ 
tions  of  motion  and  thermal  energy  until  conver¬ 
gence  is  achieved  for  each  time  step  of  the  calcu¬ 
lation. 

At  each  time  step,  the  drag,  lift,  and  moment  co¬ 
efficients  and  Nusselt  number  are  evaluated.  The 
entire  procedure  is  then  repeated  for  the  next  time- 
step.  Further  details  may  be  found  in  previous  publi¬ 
cations  of  this  research  group6-8  .  All  the  executions 
were  pursued  on  a  Dec-alpha,  a  Convex  240  and  a 
Convex  3840. 


dudvdwdT 

p  =  0’d(  =  e(‘d(  =  d(=0  it{  =  w“‘nd 


1  <  77  <  N2mid  (downstream) 


The  imposed  initial  conditions  inside  the  liquid 
droplet  are  a  quiescent  liquid  phase  and  a  uniform 
temperature  T/o  <  Tg q. 


Symmetry  Conditions 


Since  the  cylindrical  vortex  tube  advects  with  its 
axis  of  symmetry  parallel  to  the  y-axis,  symmetry 
is  maintained  such  that  we  solve  for  half  the  spheri¬ 
cal  domain  rather  than  the  entire  domain,  and  thus 
reduce  the  computational  time. 


dp  du 

dC  =  d( 


dw 

~dC 


atC  =  1,^3 


Numerical  Solution 


A  three-dimensional  implicit  finite-difference  al¬ 
gorithm  solves  the  set  of  discretized  partial  differ¬ 
ential  equations.  The  control  volume  formulation 
is  used  to  develop  the  finite-difference  equations. 
The  method  of  solution  employs  an  Alternating- 
Direction-Predictor-Corrector  (ADPC)  scheme  to 
solve  the  time-dependent  equations. 

The  overall  solution  procedure  is  based  on  a  cyclic 
series  of  estimate-and-correct  operations.  At  each 
time-step,  we  first  regard  the  solution  in  the  gas- 
phase;  the  velocity  components  are  first  calculated 
from  the  momentum  equations  using  the  ADPC 
method,  where  the  pressure  field  at  the  previous 
time  step  is  employed.  This  estimate  improves  as 
the  overall  iteration  continues.  The  pressure  is  cal¬ 
culated  from  the  pressure  correction  equation  using 
the  successive  over  relaxation  method.  The  new  es¬ 
timates  of  pressure  and  velocities  are  then  obtained. 
These  known  quantities  are  used  in  the  energy  equa¬ 
tion  to  solve  for  the  gas-phase  temperature  field. 


The  Vortex  Tube  Features 

The  vortex  is  introduced  upstream  of  the  droplet, 
advects  with  the  superimposed  uniform  flow,  and 
has  a  relatively  simple  configuration — it  is  an  ini¬ 
tially  cylindrical  tube  whose  axis  of  symmetry  is  ini¬ 
tially  normal  to  the  uniform  flow  and  parallel  to  the 
y—  axis.  The  vortex  tube  has  a  relatively  small  cen¬ 
tral  core.  Within  this  core,  the  initial  velocity  dis¬ 
tribution  in  the  vortex  tube  is  that  of  solid  body  ro¬ 
tation  reaching  an  imposed  vmax  at  radius  a.  vmax 
and  a  are  specified  at  time  t  =  0.  Outside  this  inner 
core,  the  vortex  induces  a  velocity  field  of  a  potential 
vortex;  thus,  the  velocity  induced  by  the  vortex  van¬ 
ishes  as  r  — >  00.  This  two-dimensional  vortical  tube 
is  known  as  Rankine  vortex9  ,  and  has  the  following 
stream  function10  : 


rpv(x,z,t  =  0)  =  ~^ln[(x  -  x0)2  +  (z  -  z0)2  +  <?%] 


(7) 

where  To  is  the  initial  non-dimensional  vortex  circu¬ 
lation  at  radius  <r0.  To  is  positive  when  the  vortex 
tube  has  counterclockwise  rotation,  and  xo  and  z 0 
denote  the  initial  location  of  the  center  of  the  vortex 
tube.  Note  the  vortex  tube  circulation  at  radius  a 
is  T  =  27 T<TVmax‘  More  fundamental  information  on 
the  vortex  tube  such  as  temporal  changes  in  its  tan¬ 
gential  velocity  and  vorticity  are  given  in  Ref.  [6]  . 

Flow  Interaction 

The  droplet  is  placed  in  a  uniform  flow  (here  also 
called  the  ‘base  flow*)  and  thus  gradually  develops  a 
standing  vortex  ring  in  its  aft  position.  Note  that,  in 
the  absence  of  the  vortex,  the  flow  remains  axisym- 
metric  with  respect  to  the  z— axis  (Fig.  1).  The 
vortex  is  introduced  10  droplet  radii  upstream  of 
the  droplet  and  advects  with  the  superimposed  uni¬ 
form  flow;  it  takes  about  10  residence  time  units  for 
the  vortex  to  arrive  at  the  vicinity  of  the  droplet; 
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there,  we  observe  vortex  stretching  in  the  cross- 
flow  direction  and  thus  a  full  unsteady  and  three- 
dimensional  interaction  occurs  between  the  vortex 
and  the  droplet.  The  dynamic  interaction  is  the 
strongest  when  the  vortex  is  initially  introduced  ‘on’ 
the  base  flow  symmetry  axis  of  Fig.  1;  here,  a  ‘head- 
on’  collision  between  the  droplet  and  the  vortex  is 
observed,  resulting  in  vortex  stretching  in  the  cross- 
flow  direction.  When  the  vortex  advects  ‘off’  the 
axis,  the  dynamic  interaction  between  the  two  is  rel¬ 
atively  weaker.  It  takes  nearly  25  Residence  time 
units  for  the  vortex  to  arrive  at  the  droplet,  interact 
with  it,  and  then  travel  sufficiently  far  downstream 
to  have  insignificant  influence.  Many  details  of  the 
interactions  have  been  reported  in  Ref.  [6].  Figs.  (2) 
show  how  the  velocity  and  thermal  boundary  layers, 
both  in  the  gas  phase  and  within  the  droplet  interior, 
could  be  affected  due  to  the  advection  of  the  vortex 
near  the  droplet.  Naturally,  and  as  observed,  a  more 
pronounced  departure  from  the  axisymmetric  case  is 
expected  in  the  gas  phase. 

The  Droplet  Convective  Heat  Transfer 

The  droplet  convective  heat  transfer,  represented 
by  its  Nusselt  number,  is  computed  through  Nu(t)  = 
2 afhf/kfg  (with  hl  and  kfg  being  the  convective 
heat  transfer  coefficient  and  the  gas  conductiv¬ 
ity)  which  after  a  standard  simplification  and  non- 
dimensionalization  yields 

=  ^.jinedOd* 

»(1  -T.) 

where  Ts  is  the  droplet  temperature  at  the  interface 
averaged  over  the  surface.  Since  the  cold  droplet  is 
injected  impulsively  in  the  hot  ambient  gas,  it  ini¬ 
tially  experiences  a  stronger  heat  transfer.  In  the 
base  axisymmetric  flow,  it  takes  nearly  5  resident 
time  units  for  the  droplet  Nusselt  number  to  reach 
a  steady  value.  However,  when  the  vortex  is  su¬ 
perimposed  on  the  base  flow,  the  Nusselt  number 
fluctuates  continuously  due  to  the  advection  of  the 
vortex  and  can  not  attain  a  steady  value.  It  is  there¬ 
fore  more  convenient  to  regard  overall  estimates  by 
considering  time-averaged  values  according  to 

-  1  f* 2 

Nu  =  - - -  /  Nu(t)  dt  (8) 

h  — Jn 

It  is  further  advantageous  to  normalize  Nu  us¬ 
ing  its  corresponding  value  in  an  axisymmetric  flow 
N uax ;  we  will  thus  consider  Nu/Nuax  in  our  results. 
Also,  in  the  above  estimates,  t\  =  2  and  t2  =  25; 
i.e.  we  disregard  the  data  for  t  e  (0,2);  this  is  the 


time  needed  for  the  initial  computational  fluctua¬ 
tions  in  the  pressure  drag  to  vanish.  Since  we  also 
estimate  the  normalizing  value  Nuax  after  this  ini¬ 
tial  data  exclusion,  the  final  quantity  Nu/Nuax  is 
barely  changed;  (sample  comparison  shows  the  effect 
of  the  initial  data  exclusion  in  Nu/Nuax  is  less  than 
.05%). 


3.  Results 

There  are  four  parameters  that  characterize  the 
quantitative  significance  of  the  vortex-droplet  inter¬ 
action:  the  vortex  initial  core  size  (<t0)5  its  initial 
maximum  tangential  velocity  ( v0max ),  the  offset  dis¬ 
tance  do  between  the  vortex  initial  position  and  the 
base  flow  symmetry  axis  (d0),  and  the  base  flow 
Reynolds  number  (Re).  We  place  the  initial  posi¬ 
tion  of  the  vortex  center  10  droplet  radii  upstream 
the  droplet,  either  ‘on’  the  base  flow  symmetry  axis 
(d0  =  0)  or  slightly  ‘off’  it  (d0  =  ±1,  ±2, ....).  A  pos¬ 
itive  or  negative  do  means  an  offset  distance  from  the 
z— axis  in  the  x,  z  symmetry  plane  in  the  positive  or 
negative  x  direction,  respectively;  this  is  shown  in 
Fig.  1  with  do  non-dimensionalized  by  the  droplet 
radius. 

We  specify  an  initial  radius  (<70)  for  the  vortex 
which  defines  the  vortex  core  within  which  vorticity 
is  uniformly  distributed;  we  pick  the  strength  of  this 
vorticity  so  that  the  maximum  velocity  at  the  core 
of  the  vortex  (vQmax)  represents  an  acceptable  fluc¬ 
tuation  from  the  uniform  flow.  This  fluctuation  is 
taken  to  be  less  than  the  free  stream  velocity.  For 
example,  to  represent  a  20%  fluctuation  in  the  base 
flow,  we  pick  v0max  =  .2.  Outside  the  inner  core,  the 
velocity  pattern  is  that  of  a  potential  vortex.  Thus, 
the  vortex  structure  and  strength  are  initially  fully 
characterized  by  the  two  parameters  <r0  and  VQmax, 
non-dimensionalized  by  the  droplet  radius  and  the 
strength  of  the  uniform  stream,  respectively.  Fig.  1 
shows  the  vortex  location  upstream  of  the  droplet. 
Since  in  the  absence  of  the  vortex  the  base  flow  re¬ 
mains  axisymmetric  at  all  times,  a  change  in  the 
orientation  of  the  vortex  circulation  only  rotates  the 
spatial  orientation  of  the  events;  so,  we  arbitrarily 
pick  the  counterclockwise  orientation  for  the  vortex 
in  all  our  simulations.  A  counterclockwise  rotation 
with  a  positive  offset  distance  is  the  mirror  image  of 
a  clockwise  rotation  with  a  negative  offset  distance 
so  that  clockwise  rotation  need  not  be  considered. 

In  order  to  investigate  the  droplet  heat  trans¬ 
fer  influenced  by  the  flow  fluctuations  due  to  the 
passage  of  the  vortex,  we  pursue  a  parameter 
study  to  determine  the  role  of  each  of  the  four 
above  characteristics  on  the  droplet  Nusselt  num¬ 
ber.  In  the  following  sections,  we  report  obser¬ 
vations  and  resulted  correlation  for  Nusselt  num- 
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ber  affected  by  these  four  parameters.  The  con¬ 
sidered  ranges  are  d0  =  0,  ±1,  ±2,  ±3,  ±4,  ±5; 

(To  =  .25,  .5, 1,2, 3,4;  vomax  —  .1,.2,.3, .4,  and 
20  <  Re  <  100. 

The  Effect  of  the  Offset  Distance 

Figures  3(a,b)  show  the  temporal  changes  in  Nus- 
selt  number  as  a  function  of  vortex  center  initial  po¬ 
sition  (do)  upstream  the  droplet.  The  droplet  Nus- 
selt  number  in  an  axisymmetric  flow  (i.e.  having 
the  same  Reynolds  number  and  without  a  vortex) 
is  also  shown  for  comparison.  The  substantial  dif¬ 
ference  in  the  temporal  response  is  apparent.  The 
droplet  Nusselt  number  increases  for  do  >  0  and  de¬ 
creases  for  d0  <  0;  a  vortex  with  a  counterclockwise 
circulation,  when  positioned  at  an  initial  do  >  0,  in¬ 
creases  the  relative  gas-droplet  velocity  in  the  vicin¬ 
ity  of  the  droplet  and  thus  increases  its  convective 
heat  transfer;  by  contrast,  one  with  do  <  0  decreases 
the  convective  heating  of  the  droplet. 

An  interesting  case  is  that  of  do  =  0  where  the 
droplet  Nusselt  number  goes  through  a  pattern  of 
increase- decrease-in  crease  depending  upon  the  vor¬ 
tex  location  while  advecting.  This  pattern  for  do  =  0 
is  seen  in  Figs.  3(a,b).  When  the  vortex  is  upstream 
of  the  droplet,  its  counterclockwise  circulation  in¬ 
creases  the  convective  effect,  and  thereby  the  droplet 
Nusselt  number;  very  near  the  droplet,  viscous  inter¬ 
actions  force  it  to  pass  ‘underneath’  the  droplet6  . 
(‘ underneath’  means  in  the  lower  half  of  the  x,  z  sym¬ 
metry  plane  in  Fig.  2.)  In  the  droplet  vicinity,  the 
vortex  reduces  the  convective  effect  decreasing  the 
droplet  Nusselt  number.  When  downstream  of  the 
droplet,  the  vortex  once  again  strengthens  the  con¬ 
vective  effect  and  increases  the  droplet  Nusselt  num¬ 
ber. 

Moreover,  when  the  vortex  is  initially  positioned 
on  the  base  flow  symmetry  axis  (do  =  0),  an  increase 
in  Nusselt  number  is  followed  by  a  decrease  and  thus 
the  summed  variations  yielding  Nu  are  small;  by 
contrast,  the  non-trivial  changes  in  Nu  occur  when 
do  ^  0.  The  measured  values  in  such  case  suggest 
that  Nu/Nuax  —  1  ~  tanh(do);  this  is  shown  in 
Fig.  3(c). 

Note  also  that,  in  spite  of  sensitivity  of  Nu(t) 
to  the  passage  of  the  vortex  (Fig.  3(a,b)), 
Nu/Nuax  —  1  is  nearly  invariant  to  the  vortex  for 
2  <|  do  |<  5  (Fig.  3(c)). 

The  Effect  of  the  Voi'tex  Tangential  Velocity 

Figures  4(a,b)  show  the  influence  of  vomax  on  tem¬ 
poral  Nusselt  number  when  the  vortex  advects  on 
(do  =  0,  Fig.  4(a))  or  off  (d0  ^  0,  Fig.  4(b))  the 
base  flow  symmetry  axis;  the  corresponding  pattern 
in  an  axisymmetric  flow  with  the  same  Reynolds 


number  has  been  also  shown  for  comparison.  Since 
Tv  a  Vfnax ,  vortices  with  larger  vmax  induce  a 
stronger  secondary  flow  in  the  uniform  stream;  thus, 
droplet  Nusselt  number  shows  sensitivity  to  the  vor¬ 
tex  Vo  max  • 

However,  changes  in  Nu  due  to  vomax  appear  to 
be  small  for  as  long  as  do  =  0;  (see  Fig.  4(a)).  This 
is  the  same  as  the  previous  observation  that,  as  long 
as  d0  =  0,  the  vortex  effect  on  Nu(t)  is  sometimes 
augmenting  and  sometimes  decreasing,  yielding  triv¬ 
ial  changes  from  the  axisymmetric  value  in  Nu.  This 
was  also  seen  when  studying  the  effect  of  do.  The  ob¬ 
served  negligible  changes  in  Nu  should  not  be  inter¬ 
preted  as  the  insensitivity  of  the  droplet  heat  trans¬ 
fer  to  the  passage  of  the  vortex  or  the  vortex  vomax ; 
simply  stated,  at  do  =  0,  we  find  cancellations  in  the 
averaging  process.  This  cancellation  is  not  encoun¬ 
tered  in  computing  Nu  when  do  ^  0,  such  as  the 
cases  shown  in  Fig.  4(b);  later,  we  will  discuss  this 
more  interesting  effect  of  vortex  velocity  that  occur 
when  do  ^  0. 

Data  indicate  that,  overall,  Nu/Nuax  —  1  nearly 
follows  a  linear  dependence  on  vomax. 

The  Effect  of  the  Vortex  Radius 

Figures  5(a,b)  show  the  influence  of  (Tq  on  tem¬ 
poral  Nusselt  number  when  the  vortex  advects  on 
(d0  —  0,  Fig.  5(a))  or  off  (do  /  0,  Fig.  5(b))  the  base 
flow  symmetry  axis;  the  corresponding  pattern  in  an 
axisymmetric  flow  with  the  same  Reynolds  number 
has  been  also  shown  for  comparison. 

Analogous  to  the  previous  observation,  since  Tv  oc 

vortices  with  larger  radii  introduce  stronger  sec¬ 
ondary  flow  in  an  otherwise  uniform  stream  near  the 
droplet;  thus,  droplet  Nusselt  number  appears  sensi¬ 
tive  to  change  in  the  vortex  radius.  When  the  vortex 
advects  on  the  base  flow  symmetry  axis  (do  =  0), 
in  spite  of  the  clear  temporal  sensitivity,  the  time- 
average  values  of  Nusselt  number  compared  to  those 
in  an  axisymmetric  flow  (having  the  same  Reynolds 
number)  change  by  only  .5%  when  v0max  =  2  and 
by  about  1%  when  VQmax  =  A  so  that  we  may  note 
Nu/Nuax  —  1  ±  1%.  Similar  to  the  case  in  the  pre¬ 
vious  section,  the  trivial  net  change  in  the  computed 
Nu  in  such  cases  is  due  to  a  combined  effect  of  do  =  0 
in  the  simulation  and  the  nature  of  time-averaging 
of  Eqn.  (8);  it  should  not  be  interpreted  as  the  in¬ 
sensitivity  of  the  droplet  convective  heat  transfer  to 
the  passage  of  the  vortex,  (see  Fig.  5(a)). 

Thus,  when  studying  the  effect  of  (Tq  on  Nu/Nuax , 
larger  values  are  resulted  when  the  vortex  advects  off 
the  base  flow  symmetry  axis,  (Fig.  5(b));  later,  we 
shall  discuss  this  effect  of  vortex  radius  at  do  ^  0. 
Data  indicate  that  Nu/Nuax  —  1  has  a  slightly 
stronger  dependence  on  <r0  than  it  does  on  Vomai. 
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The  Effect  of  the  Flow  Reynolds  Number 

Figures  6(a,b)  show  the  influence  of  the  base  flow 
Reynolds  number  on  the  temporal  Nusselt  number 
when  the  vortex  advects  on  (do  =  0,  Fig.  6(a))  or  off 
(do  i1  0,  Fig.  6(b))  the  base  flow  symmetry  axis;  the 
corresponding  pattern  in  an  axisymmetric  flow  with 
the  same  Reynolds  number  has  been  also  shown  for 
comparison. 

The  higher  the  Reynolds  number,  the  stronger 
the  heat  transfer  fluctuations  induced  by  the  vor¬ 
tex,  even  though  the  strength  of  the  vortex  re¬ 
mains  the  same;  thus,  for  example,  a  30%  fluctu¬ 
ations  (i.e.  vomax  =  -3)  in  a  uniform  flow  with 
Re  =  100  has  a  stronger  effect  on  the  droplet  Nus¬ 
selt  number  than  the  same  fluctuation  does  in  a 
flow  with  Re  —  20.  The  reason  is  simple:  flows 
with  lower  Reynolds  numbers  are  relatively  more 
viscosity-dominated  and  thus  vortex-induced  iner¬ 
tial  changes  are  relatively  more  damped. 

Analogous  to  the  previous  observations,  when 
do  =  0  (Fig.  6(a)),  there  is  barely  change  in  time- 
averaged  values  of  Nusselt  number  due  to  change  in 
base  flow  Reynolds  number;  Nu/Nuax  —  1  changes 
over  the  range  of  Reynolds  number  by  only  .1%  when 
vomax  =  .2  and  by  .5%  when  v0max  =  .4.  Again,  this 
is  mostly  the  consequence  of  the  choice  do  =  0  and 
the  nature  of  time-averaging  in  Eqn.  (8).  Naturally, 
when  dp  ^  0  (Fig.  6(b)),  the  computed  changes  in 
Nu/ Nuax  —  1  are  stronger;  such  effects  are  seen  in 
Fig.  6(b)  and  will  be  discussed  in  the  next  section. 

Data  indicate  that  Nu/Nuax  -  1  ~  Re  0  note 
since  Nuax  ~  Rea ,  where  the  exponent  ranges  from 
•427  to  .573  (see  Appendix),  so  that  approximately 
Nu  —  Nuax  ~  Re.  This  implies  that  the  time- 
averaged  perturbations  in  Nusselt  number  have  a 
stronger  dependence  on  Reynolds  number  than  Nus¬ 
selt  number  itself  does. 
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within  the  range  of  our  parameter  study,  20  <  Re  < 
100;  .25  <  (70  <  4,  .1  <  t>0maar  <  .4,-5  <  do  <  5. 
All  our  simulations  are  for  an  n-octane  droplet  and 
so  Pr\  =  8.527.  (For  all  values  of  Re  in  the  given 
range,  we  have  verified  the  above  correlation  for 
I  do/<?o  |<  2;  work  is  under  way  to  confirm  fully 
the  correlation  for  |  d0/cr0  |>  2,  as  well.)  Modest 
modification  (  <  10%)  to  the  exponents  is  made 
by  a  least-squares  fit.  The  correlation  coefficient 
for  the  above  fit  is  .985.  A  plot  of  this  correla¬ 
tion  and  its  comparison  with  our  data  is  shown  in 
Figs.  7(a,b). 

The  significance  of  the  our  observation  should  not 
be  underestimated  for  the  following  reason:  note 
that  the  problem  at  hand  is  unsteady  by  nature; 
by  contrast,  it  has  been  well-established  that  ma¬ 
jority  (but  not  all)  of  self-similar  patterns  in  fluid 
flows  and  their  coupled  heat  transfer  mechanisms  es¬ 
tablish  themselves  when  the  transient  response  has 
vanished  and  only  when  the  equilibrium  stage  has 
established  itself.  The  observation  made  from  Eqn. 
(10)  and  Figs.  7(a,b)  reveals  however  that,  in  spite 
of  the  unsteady  nature  of  vortex- droplet  interaction, 
self-similar  correspondence  nevertheless  exist  in  the 
time-averaged  behavior. 


Global  Self- Similarity 

In  the  absence  of  the  vortex,  Nu  and  Nuax  are 
identical,  and  one  may  approximate  the  Nusselt 
number  correlation  from  Ref.[  11]  for  a  spherical 
droplet  in  an  axisymmetric  flow  by  using  the  solid 
sphere  correlation4 

Nu  =  1  +  (1  +  Pr  Re)1*3 Re-077  (9) 

or  one  may  use  the  new  correlation  for  axisymmetric 
flow  past  liquid  spheres  presented  in  the  Appendix. 
With  the  vortex  present  near  the  droplet,  however, 
the  axisymmetric  flow  correlation  loses  its  applica¬ 
bility  and  a  new  correlation  accounting  for  the  effect 
of  the  advecting  vortex  on  the  droplet  heat  transfer 
is  to  be  used. 

We  have  produced  a  correlation  that  includes  such 
effects;  all  our  data  collapse  into  the  functional  form 


4.  Conclusions 

We  have  investigated  the  unsteady  interaction  be¬ 
tween  a  cylindrical  vortex  tube  and  a  liquid  droplet 
in  a  uniform  flow  by  pursuing  a  study  of  each  of 
the  four  parameters  affecting  the  strength  of  this 
interaction  and  the  modification  of  the  droplet  Nus¬ 
selt  number.  We  have  paid  particular  attention  to 
the  effect  of  these  parameters  on  the  transient  and 
time-averaged  values  of  the  droplet  Nusselt  number. 
A  correlation  quantifying  the  effect  of  the  advecting 
vortex  on  the  droplet  heating  has  been  produced. 
The  reported  correlation  compliments  the  existing 
ones  for  droplet  heating  in  axisymmetric  flows  that 
occur  in  the  absence  of  an  advecting  vortex  11  . 

When  the  vortex  advects  towards  and  then  past 
the  droplet  starting  upstream  ‘on’  the  symmetry 
axis  of  the  base  flow,  the  droplet  Nusselt  number 
first  increases  and  then  decreases  in  time;  thus,  the 
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time- averaged  Nusselt  number  is  nearly  equivalent 
to  that  in  an  axisymmetric  flow  even  when  the  vor¬ 
tex  simulates  up  to  40%  fluctuation  in  the  base  flow. 
Conversely,  when  the  vortex  advects  ‘off’  the  base 
flow  symmetry  axis,  time-averaged  Nusselt  number 
is  changed  due  to  the  vortex- induced  fluctuations; 
whether  this  has  an  augmenting  or  inhibiting  effect 
on  the  droplet  heating  depends  on  two  factors:  the 
vortex  circulation  orientation  and  also  whether  it  ad¬ 
vects  ‘above’  or  ‘below’  the  symmetry  axis  in  the 
plane  of  symmetry.  Time- averaged  Nusselt  num¬ 
ber  monotonically  varies  with  the  ratio  of  the  vortex 
initial  offset  distance  to  the  vortex  initial  core  size 
varying  from  0  to  2  (i.e.  0  <|  (Iq/<tq  |<  2);  when 
this  ratio  is  greater  than  2,  time-averaged  Nusselt 
number  remains  nearly  unchanged  (Figs.  7). 

In  the  range  0  <|  cfo/co  |<  2,  change  in  the  droplet 
Nusselt  number  has  a  linear  dependence  on  the  vor¬ 
tex  circulation,  a  weak  dependence  on  the  vortex 
initial  core  size,  and  ~  tanh(do/<ro).  We  believe  this 
occurs  because,  when  the  vortex  arrives  at  the  vicin¬ 
ity  of  the  droplet,  the  droplet  is  encompassed,  partly 
or  wholly,  within  the  vortex  inner  core;  thus,  the  per¬ 
turbations  in  flow  configuration  depend  on  both  the 
vortex  core  size  and  its  circulation.  (Dependence 
on  the  vortex  initial  maximum  tangential  velocity 
is  embedded  within  the  dependence  on  the  vortex 
circulation.) 

The  situation  is  different  when  |  do/cr0  |>  2  :  when 
the  vortex  arrives  at  the  vicinity  of  the  droplet,  the 
droplet  is  totally  outside  the  vortex  inner  core  and 
thus  it  is  ‘only’  the  vortex  circulation  that  matters 
to  the  droplet.  Thus,  it  is  not  surprising  that  the 
computed  fluctuations  in  this  range  in  Eqn.  (10) 
only  depend  on  the  vortex  circulation. 

We  emphasize  that,  in  parallel  studies  investi¬ 
gating  changes  in  fluid  dynamic  properties  of  solid 
spheres  (lift  and  moment  coefficients)6  ,  similar  de¬ 
pendence  of  fluctuations  in  fluid  dynamic  properties 
on  vortex  parameters  have  been  observed.  When  the 
vortex  advects  near  the  sphere,  the  induced  fluctua¬ 
tions  appear  to  depend  both  on  the  vortex  size  and 
its  circulation;  by  contrast,  with  the  vortex  advect- 
ing  relatively  far  from  the  sphere,  the  changes  ap¬ 
pear  to  depend  only  on  vortex  circulation.  Marked 
dependence  of  Nusselt  number  fluctuations  on  the 
fluid  dynamics  variations  explains  the  similarity  be¬ 
tween  the  conclusions  of  the  two  studies. 

A  rather  surprising  observation  here  is  that, 
while  in  the  absence  of  vortical  structures  Nuax  & 
0(Re  0  5)  thus  prompting  one  to  anticipate  its  per¬ 
turbations  Nu  -  Nuax  «  0(Re  °  5)  as  well,  our  ob¬ 
servations  here  show  that  this  is  not  the  case;  in¬ 
stead,  the  induced  perturbations  follow  «  O(Re)  . 

Based  on  our  findings,  it  is  speculated  that,  in  a 
spray  combustion  system,  the  vortex-droplet  inter¬ 
actions  within  the  Kolmogorov  scale  can  have  signif¬ 


icant  effects  on  the  droplet  convective  heat  transfer. 
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Appendix 

In  an  axisymmetric  flow,  and  following  the  as¬ 
sumption  that  the  effect  of  the  droplet  internal  cir¬ 
culation  on  the  droplet  heating  is  small,  one  may  use 
the  Nusselt  number  correlation  from  Ref.[  11]  for  a 
rigid  sphere 

Nu  =  1  +  (1  +  Pr  Re)1/3 Re0  077  (11) 

for  Re  <  400.  This  correlation  has  been  recom¬ 
mended  for  spray  combustion  computations4  .  For 
particles  having  Re  >  10  in  media  with  Pr  &  1,  this 
suggests  Nu  —  1  ~  Re0  410. 

We  have  modified  this  suggested  correlation  by 
accounting  for  the  droplet  internal  circulation.  Our 
calculations  (n-octane  droplet  in  air,  10  <  Re  <  100) 
suggest  that  a  fit  of  the  form 

Nu  —  1  =  .927  Re  0  427  10  <  Re  <  100.  (12) 

The  authors  note,  however,  that  this  form  is  some¬ 
what  clumsy  since  one  should  expect  Nu  =  2  at 
Re  =  0.  A  preferred  fit  therefore  is 

Nu  -  2  =  .412  Re  0  573  10  <  Re  <  100.  (13) 

Each  of  these  fits  has  a  correlation  coefficient 
greater  than  .99. 

It  is  suggested  that  the  slight  increase  in  the  expo¬ 
nent  of  Re  (compared  to  the  exponent  .410  resulted 
from  Eqn.  (11))  is  a  contribution  of  the  droplet  in¬ 
ternal  circulation.  That  is,  the  boundary  layer  and 
thermal  layer  thicknesses  are  decreased  slightly  due 
to  the  motion  along  the  interface.  This  results  in  an 
increase  in  heat  transfer  rate  from  the  hot  ambient 
gas  to  the  cold  droplet. 
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Fig.  1.  Flow  geometry  and  coordinates. 


Fig.  2.  (Previous  page.)  Influence  of  an  advecting  vortex  on  the 
layers  of  a  moving  cold  droplet  in  a  hot  gas  (d0  =  2,  a0  =  4,  v0 


(a)  gas  phase  velocity  field,  axisymm.  case. 

(b)  same,  affected  by  an  advecting  vortex. 

(c)  gas  phase  thermal  field,  axisymm.  case. 

(d)  same,  affected  by  an  advecting  vortex. 


velocity  and  thermal  boundary 
*  =  A,  Re  =  100,  t  =  10.5). 

(e)  liquid  interior  velocity  field,  axisymm.  case. 

(f)  same,  affected  by  an  advecting  vortex. 

(g)  liquid  interior  thermal  field,  axisymm.  case. 

(h)  same,  affected  by  an  advecting  vortex. 
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Figs.  3.  (Previous  page.)  Influence  of  the  vortex  initial  offset  distance  ( d0 )  on  Nu(t): 

(a)  vortex  advecting  with  d0  in  —  5  <  d0  <  5,  (a0  =  1  ,Re  =  100,  v0max  =  .2); 

(b)  same,  Vomai  =  _ 

(c)  appearance  of  Nu/Nuax  —  1  ~  tanh(d0)- 
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Figs.  4.  Influence  of  the  vortex  initial  tangential  velocity  (v0max)  on  Nu(t)  : 

(a)  vortex  advecting  with  d0  =  0,(iie  =  100,  a0  =  1):  ; 

(b)  same,  d0  =  1. 
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Figs.  5.  Effect  of  the  vortex  initial  radius  («r0)  on  Nu(t )  : 

(a)  vortex  advecting  with  d0  =  0,  (Re  =  100,uOmax  =  -2); 

(b)  same,  d0  =  1. 


